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Dynamics of a dissociating gas 
Part I Equilibrium flow 


By M. J. LIGHTHILL 
Department of Mathematics, University of Manchester 


(Received 12 October 1956) 
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1. INTRODUCTION 

Atmospheric dissociation will be appreciable in the neighbourhood of 
projectiles travelling at speeds greater than 2 km/sec. ‘This introductory 
paper on possible effects of dissociation on the airflow, and hence on the 
drag, stability and aerodynamic heating of such projectiles, is intended 
mainly as a source of ideas for later, more comprehensive investigations. 

The problem of incorporating the effects of the large energy change 
involved in dissociation into the standard theory of gas flow appears at the 
same time so important and so formidable that it is worth approaching 
slowly. One may usefully begin, on both the theoretical and experimental 
sides, by eliminating the less essential complications which arise from the 
detailed composition of air, and studying the dynamics of a pure dissociating 
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diatomic gas. By keeping within bounds in this way the difficulty of 
imagining what is happening in the flow, one reduces the risk that some 
important feature will be missed. 

The theoretical treatment of the dynamics of a pure dissociating gas 
which follows is concerned principally with properties at densities between 
10° and 1 of N.T.P., and at temperatures such that dissociation occurs 
appreciably but ionization is negligible. For the gases, O, and N,, which 
will be principally used for numerical illustration, this means temperatures 
in the approximate range 3000 K to 7000°K. Other gases suitable for 
pioneer experimental work on the problem and comparison with theory are 
the halogens, which already dissociate appreciably at 1000° K and so can 
be studied with lighter equipment. ‘Thus, Britton, Davidson & Schott 
(1954) have used I,, and Palmer (1955) Brg. 

The problem falls into three main parts, as follows. 


(i) ‘The theory of flow in which thermodynamic equilibrium is main- 
tained everywhere ; this is a good approximation to flow outside the boundary 
layer if the time scale for flow past the body is large compared with the time 
scale tor dissociation or recombination. ‘This theory is given in § 3, after 
the necessary equilibrium thermodynamics has been discussed in § 2. 

(ii) ‘he quasi-equilibrium theory of the transport properties, including 
radiative heat transfer as well as convective heat and momentum transfer 
in the boundary layer; this theory (to be given in Part II) is a good 
approximation if the time scale for diffusion across the boundary layer is 
large compared with the time scale for dissociation or recombination. 
Actually, the time scales for flow past the body and for diffusion across 
the boundary layer are the same (which is, indeed, what determines the 
thickness of the boundary layer), and so there is just one condition for 
theories (1) and (ii) to be valid. 

(iii) ‘he theory of flow in which this condition is not satisfied, so that 
large departures from equilibrium occur (not just the small ones supposed 
under heading (ii) to be responsible for the heat and momentum transport) ; 
this theory (to be given in Part III) includes the theory of the extended 
character of the shock wave and the effect of this on the flow behind it, 
and also the effect on boundary-layer behaviour of the delay in recombination 
of free atoms diffusing into the neighbourhood of the cooled wall. 

The paper begins with a study (§ 2) of the equilibrium statistical thermo- 
dynamics of a pure dissociating gas. Here, an approximation is found 
which greatly simplifies the analysis and yet introduces only small errors 
for particular gases. One may speak of a hypothetical gas satisfying this 
approximate form of the thermodynamic equations as an ‘ideal dissociating 
gas’. ‘The properties of an ideal dissociating gas are given completely once 
three constants 7’, p, and u, (the ‘characteristic’ temperature, density and 
specitic energy for dissociation) are given. Different pure dissociating gases 
differ in the values to be assigned to these constants, but otherwise (to the 
good approximation with which they can be treated as ‘ideal’ dissociating 
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gases) have all the same thermodynamics; a single theory, therefore, will 
do for all, whereas if the aproximation were not made separate calculations 
would need to be done for each gas. 

Approximations equivalent to those in the thermodynamical theory of 
an ideal dissociating gas are made also in the discussions of quasi-equilibrium 
flow and flow involving large departures from equilibrium in Parts II and III. 
Physically, they amount always to taking the vibrational modes of motion 
of the molecules as everywhere excited to just half their full ‘classical’ 
energy content. At the high temperatures at which this would be expected 
to be a serious underestimate, the molecules have been so reduced in number 
by dissociation, and the energy absorbed in this process has been so large 
a fraction of the total internal energy, that the latter is in fact underestimated 
by a very few percent. Under these circumstances the gain in simplicity well 
compensates for the loss in accuracy, as has been found many times before 
when hypothetical ideal fluids have been introduced in hydrodynamics. 

Another reason why in this paper the author has avoided the complications 
due to the variation in degree of excitation of the vibrational mode of motion 
is that these complications were rather fully treated (in the absence of 
dissociation) in his recent survey article in the G. I. ‘Taylor Anniversary 
Volume (Lighthill 1956). In the present paper the enquiry into the effects 
of molecular constitution on gas dynamics, begun in that article, is continued 
with a study of the effects of dissociation, but to simplify matters the 
concomitant effects of variation of vibrational excitation are here excluded. 

The papers and books listed in the bibliographies have greatly assisted 
the author in preparing the three parts of this paper. However, like his 
other papers, it could never have been written without extensive private 
discussion and correspondence with friends. ‘The help of Mr E. Wild 
has again been invaluable, and it is a pleasure to thank in addition 
Dr W. Chester, Mr N. C. Freeman, Dr A. G. Gaydon, Dr W. C. Griffith, 
Dr A. R. Kantrowitz, Dr L. Lees, Mr D. J. Lyons, Dr J. S. Rowlinson, 
Dr k. Stewartson, Dr P. Thompson, Mr A. K. Weaver, Dr M. D. Van Dyke, 
and Mr H. K. Zienkiewicz for many useful suggestions. 


2, EQUILIBRIUM THEORY OF A PURE DISSOCIATING DIATOMIC GAS 
2.1. Law of mass action for dissociative equilibrium 
The law governing the equilibrium concentrations in a dissociative 
process A,=A+A (1) 
is well known (see for example Fowler & Guggenheim 1939, §§ 502-508). 
We have 
2) 
where m, is the number of A atoms in volume V at a given temperature 7, 
and f, is the partition function of A, namely the sum 
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over all quantum states of the A atom in the volume V; n,, and f, are 
defined similarly, but the energies « in both f, and f,, are commonly 
measured from the energy of the atom and molecule respectively when 
at rest in their ground states. ‘The factor e?/*” (where D is the dissociation 
energy per molecule) is therefore attached to f, to bring both sets of 
energies to a common origin. Equation (2) says really that the concen- 
trations will adjust themselves in direct proportion to the number of states 
available at the given volume and temperature, bearing in mind the reduced 
availability of high-energy states which results from the Maxwell—Boltzmann 
distribution law. 

The theory as described neglects complications due to gas imperfection, 
that is, contributions to the energy at any given instant from such interactions 
between molecules as are happening at that instant, which are associated 
with pressure reduction due to attraction between molecules and pressure 
increase due to overcrowding (that is, to reduced availability of translational 
states). The corrections for these effects, whose expression by means of 
a second virial coefficient may be extrapolated above the temperatures at 
which measurements have been made by the use of a Lennard-Jones (12, 6) 
potential (see for example Hirschfelder, Curtiss & Bird 1954) are of 
magnitude about 1/p or less, if the density p is expressed in gm/cm’. 
We shall not be concerned in this paper with densities at which such a 
correction is appreciable. 

The quantum states to be enumerated in (3) are combinations of 
translational, rotational, vibrational and electronic states, which can be 
regarded for practical purposes as independent of one another. (Nuclear 
states may be ignored, as any contribution to the partition function due 
to them will cancel in the ratio (2).) Using indices 7, R, V and E for the 
four kinds of state, we have therefore 


(4) 


2.2. Detailed expressions for partition functions 

The translational partition function f7 is found by replacing the sum (3) 
by an integral over the classical ‘ phase space’, cells of which of volume h° 
contribute one state each to the sum, as 


where m is the mass of the 4 atom; clearly, f, is the same with m replaced 
by 2m. ‘The rotational partition function is similarly found (a factor } 
being inserted, however, to allow for the indistinguishability of the atoms) 
to be 


(6) 


= AT) 


where / is the moment of inertia of the 4, molecule. Equation (6) defines 
a ‘characteristic temperature for rotation’, 7. The summand in (3) varies 
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gradually enough for the sum to be replaced by an integral as in (6) provided 
that T >T,. For O,, T, = 2-07°K, and for N,, 7, = 2-86° K, so that 
this requirement is amply satisfied in the present application. 

The lower vibrational states of A, differ in energy from the gound state 
by 0, hv, 2hy, ..., where v is the frequency calculated classically from the 
mass m of the A atom and the curvature of the potential-energy curve at 
equilibrium. The vibrational partition function computed by assuming 
that the arithmetic progression of energy levels continues unchanged to 
infinity is 

1 1 
fh, ~ ~ (7) 


where 7. = hv/k, the characteristic temperature for vibration, is 2230°K 
for O, and 3340°K for N,. More accurately, one may calculate f” from 
the Morse (1929) energy levels of a molecule with dissociation energy D 


f= > exp| 1- (3) 


O0<n<2Dihv 


as 


but the ditference from (7) can be shown to be small when 7’ < D/k, which 
is 59000° K for O, and 113000° K for N,. (The same condition is amply 
adequate for the absence of any significant interaction between the rotational 
and vibrational degrees of freedom.) 

The electronic partition functions f{ and f4 consist normally of only a 
very few terms, since the higher electronic states are filled to an extent 
negligible from the thermodynamic point of view (although not from the 
spectroscopic). For oxygen, we may take 

if terms of order e?°°"" can be neglected. ‘The values quoted arise from 
the fact that the ground state of oxygen is a triplet, states of weight 5, 3 
and 1 (that is, spectroscopic terms *P,, ?P, and 3P)) being energetically 
very close to one another; while the ground state (#2) of O, is of weight 3 
and the next state (1A), of weight 2, is not quite negligibly far away. (Note 
that at the temperatures (> 2000°K) at which dissociation of O, 
appreciable f% is practically 9.) 

For nitrogen, the electronic partition functions are still simpler, and 

= 4, FR, = 1, (10) 
if terms of order e~*8°? can be neglected. 

Note that all the characteristic temperatures which have been quoted 
are known with considerable accuracy from absorption and emission spectra 
associated with transitions between the various states of the atom or molecule. 
Numerical values used here are all taken from Fowler & Guggenheim (1939) 
and Gaydon (1953). The only value which has been controversial until 
recently is the dissociation energy of N,, but the value used (9-76eV per 
molecule) is now firmly established. 
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4 2.3. Equations of equilibrium in a hydrodynamically convenient form 
For hydrodynamical calculations it is convenient to express the results 
in terms of the density, namely 


m(n,+2n,,) 


(11) 
and the proportion by mass of A atoms in the mixture, namely 
ns 
4= — 12 


Then equation (2), with all the results of equations (4) to (10) substituted 
in it, becomes 


a? n*, m 
2n,(n,+2n,)  2pV 


n 
| _ eT yT (f4)° 


he 


— Pd (13) 

p 
where 7\, = D/k is a characteristic temperature for dissociation, and p,, 
defined as the contents of the curly brackets in (13), is a characteristic 


density. 
TiC K) | o,(gm'cm"*) when K) is 
5 | | 1000 2000 3000 4000 5000 6000 7000 
—_—} 
| Oxygen | 59000 145 170 166 156 144 133 = 123 
| Nitrogen | 113000 | 113 135 136 133 128 Ws 118 


‘Table | 


‘Table 1 gives 7\,, and also p, as a function of 7, for both oxygen and 
nitrogen. It will be seen that the variation of p, over the large temperature 
range from 1000° to 7000° is very slight for both gases, especially when 
compared with the enormous variation of e-7¢7 in this range, and for 
practical purposes the useful simplification of regarding p, as a constant 
(say, 150 for oxygen and 130 for nitrogen) should lead to negligible errors. 

. This approximation will be made throughout what follows. 
i, The rather close agreement between the values of p, for the two gases 
. is due to an accidental cancelling: the abundance of rotational and 
vibrational states in O, as compared with N, would make p, substantially 
smaller for oxygen than for nitrogen, were it not that the abundance of 
electronic (and, to a minor extent, of translational) states in O as compared 
ee with N works the other way. 
Now, for atmospheric values of the density p, p,/p is at least 10°. This 
consequence of the far greater number of translational states available to 
the gas as A than as A, at these densities explains why dissociation first 
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sets in at temperatures rather small compared with the characteristic 
temperature 7',. Thus, for p,/p = 10°, (13) shows that « is already 0-05 
(5°, dissociation) when 7/7, = 0-057, and « is 0-95 (95°, dissociation) 
when 7/7, = 0-116. For densities typical of the upper atmosphere, with 
(say) p,/p = 10°, these values of 7/7, would be reduced to 0-045 and 0-076 
respectively. Figure 1 is a diagram in which by use of a reciprocal tem- 
perature scale it has been made easy to read off the degree of dissociation 
« for any values of 7/7, and p/p,. 


9, 
(3) = 7 6 5 


T/T, 
0-035 0.0375 0:04 0-055 0:06 0065007 0:08 obs O10 O12 O14016 0-20 


Figure 1. Plain lines give « (the proportion by mass of free atoms in an ideal 
dissociating gas) against 7’7T, for log,(p,/p) = 7,6,5. The diagram is 
such that linear interpolation with respect to logi(p,/p) at any horizontal 
level is accurate. For oxygen, p, — 150 gm/cm* and T, = 59000° K. For 
nitrogen, p, = 130gm/cm* and 7, ~ 106000° K. Broken lines are lines 
of constant i/u 4, the enthalpy divided by the heat of dissociation. For these, 
linear interpolation in the vertical direction is almost accurate. 


We shall need also the equation of state of the gas mixture. The pressure p 

is equal to the volume density of translational momentum transfer in a 

given direction per unit time, which occurs at an average rate RT per atom 

or molecule, giving 

V om? ( +). ( ) 

Here k/2m is often written as R; it is the gas constant per gram of A). 

(Multiplied by 1+., it is the gas constant per gram of A and A, mixed 

in the proportions «:(1—«) by mass.) It may be more convenient here 
to introduce a characteristic pressure for dissociation, 


k D 
Pa = (=) Pa Ty = (15) 
in terms of which (14) becomes 
T 
(£)(7-)a+». (16) 


Values of p, are given in table 3 below. 


Cr \ \ 
x \ ua 
a \ 
\ 2-0 
\ 
1-25 
L ~ 
0-4 
~N 
| 0-75 
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We need also an expression for the internal energy as a function of 
temperature and composition. ‘This is obtained from the partition functions 
of §2.2. ‘The total internal energy U in volume V is 


) 
U = ng sp log fa, (17) 
since the average energy per A atom is 
= k7 ap log tar (18) 


and similarly with A,; the term } Dn, is needed because in f, and f,, the 
zeros of energy are different, two atoms of A at their zero of energy having 
energy D more than a molecule of A, at its zero. 

Now, to be consistent with the approximation p, = constant suggested 
above (which implies f3, « f4,) the curly bracket in (17) should be taken 
proportional ton, +2n,, and since f, is almost exactly proportional to T??, 
the coefficient of n , + 2n,, might reasonably be taken as 3/(27), leading to 
the approximation 

U = 3kT(n,+2n,,)+3Dn,. (19) 
‘This could be physically interpreted by saying that atoms and molecules 
both have average translational energy $k7, but that molecules have in 
addition an average rotational and vibrational energy of about 3kT (it 
varies at most between k7 and 2k7), electronic contributions being neglected. 
Table 2 gives the accurate values of the coefficients of n, and 2n,, in 
(U—!Dn,)/kT, obtained from (17), for both oxygen and nitrogen, for 
comparison with 1-5, the approximate value of each coefficient suggested 
by (19), At first sight the variations might be thought important, but in 


| 
T( BK) | 1000 2000 3000 4000 5000 6000 7000 
Coefficient Oxygen | 160 1:55 1°54 1:53 1-52 1°51 
of n4 Nitrogen 150 1:50 1°50 1:50 1-50 1°50 1-50 
Coefficient Oxygen | 1:38 1°54 1:62 1:68 1:75 1:77 
of | Nitrogen | 1:31) 1:44 1:52 1:57 1:60 1:62 1-64 
Table 2 


practice they are swamped by the term involving the large dissociation 
energy D. Accordingly (19) will be used throughout what follows. It 
gives for the energy per unit mass u, commonly used in hydrodynamics, 
the expression 


He 


U 3k D 
20 
In terms of a characteristic energy per unit mass 


D 
om (21) 


u 


| 


of 
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(the dissociation energy per gram of A,), equation (20) becomes 
(22) 


ug 
Values of p,, u, and u!/? for oxygen and nitrogen are given in table 3. 


pq (atmospheres) | (k cal'gm) | (cm*’sec*) | v, = u}/? (km/sec) 


Oxygen 10" 3-67 3-9 
Nitrogen 4:1 «10° °8-02 335 x 10> 5:8 
Table 3 


Here, uj’? (= v,, say) is a characteristic velocity for dissociation; fluid 
with this velocity has enough kinetic energy to provide half the energy 
required to dissociate it completely. Note also that p, = pgv? so that p, 
is typical of the pressures obtained by stopping a flow at velocity v, and 
density pz. It exceeds the pressures which will be found in a real flow 
with velocity z,, by the same large factor by which p, exceeds the real 
density p. 


2.4. Summary of thermodynamics of an ideal dissociating gas 


We shall use the expression ‘ideal dissociating gas’ to denote a gas for 
which equations (13), (16) and (22) hold with p, (and hence also p,) replaced 
by a constant. For calculations concerning changes in such a gas between 
different states of equilibrium, it may be convenient to use 7°), p,, p, and uz 
as the units, respectively, of temperature, density, pressure and specific 
internal energy (in a flow problem this implies using u)/°= v, as the unit 
of velocity). In these units the fundamental equations (13), (16) and (22) 
take the simple forms 

23 
» u=3T+a. (23) 

The specific enthalpy, which plays an important part in steady flow 

problems, has the value 


when measured with u, as unit. Curves of constant 7 are included in 
figure 1, so that the enthalpy required to produce a given degree of 
dissociation at a given density can be read off conveniently. 

The specific entropy § will also be needed. It is obtained most simply 
from the relation 


du+pd(p)  3dT +dx+ pT(1+2)d(p~) 


dS = 
T 
3dT (1 +a)dp- (25) 
px 
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Hence, by integration, 
S = 3log T+ «(1 —2logx)—(1—z)log(1—x)—(1+z)logp+const. (26) 


The value of the constant can be found from the detailed expressions for 
the partition functions, but will not be needed in the present paper, since 
no chemical changes other than dissociation and reassociation are considered. 
Whenever numerical values of S are given, the constant in (26) has simply 
been taken as zero. Note that the unit in which S is measured in (26) is 
u,/ 1, = k/2m, simply the gas constant per gram of A). 

The properties of the ideal dissociating gas may be conveniently 
summarized in a classical enthalpy-entropy diagram (Mollier diagram) 
such as figure 2, which shows the lines of constant temperature, pressure, 
density and composition in the region of the diagram corresponding to 
pressures between (3 « 10%) and (3 x 10~*) times the characteristic pressure. 
‘These were obtained by selecting pairs of values of « and p, determining 7 
from (23) and then 7 and S from (24) and (26), and also by selecting pairs 
of values of p, 7, determining x as (1+p7-exp 7-')-1? and using (23), 
(24) and (26) to obtain 7 and S as before. Note how the enthalpy on the 
isothermals rises as the gas dissociates with decreasing pressure. ‘This 
sloping up of the isothermals causes the isobars to be unusually straight 
in the region of dissociation, since the slope di/dS = T of an isobar varies 
only slowly along it. The lines of constant density show a similar behaviour. 


4 
20 2s 30 


Figure 2. Enthalpy-entropy diagram for ideal dissociating gas (units as in § 2.4). 


It is outside the scope of this paper to consider the behaviour of 
dissociating gases in internal aerodynamics, and therefore the ‘Fanno’ 
lines representing flow in ducts of constant area have not been included 
in figure 2, though they could easily be obtained. One may remark, however, 
that the properties exhibited in figure 2 make a dissociating gas an attractive 
possibility for heat engines. ‘Thus, a closed-cycle gas turbine operating in 
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a large temperature range, which may be desirable for utilizing the energy 


output of a nuclear reactor, may in principle reach a high thermal efficiency 
with a gas in which the temperature increases so slowly along the isobars. 
At the sort of temperatures that could be conveniently dealt with, a gas 
like I, with 7, = 18000°K might be suitable. By operating up to 0-08 
or 0-09 of this characteristic temperature one might achieve the favourable 
general shape of the steam-turbine cycle with a temperature range many 
times as great. 

Note that, as «->( at low temperatures, the ideal dissociating gas 
becomes a perfect gas with constant specific heats, in the ratio y = c,/¢, = 4. 
(Physically, this is because the vibrational degrees of freedom are taken to 
be half-excited even at low temperatures.) It might be thought desirable 
for the ideal dissociating gas to have y » = at low temperatures, but, since 
projectiles travelling at speeds like those mentioned in § 1 cannot in practice 
have their surface kept below 1000°K, 4 is a reasonably realistic low- 
temperature limit for y in the field of flow around such a projectile. 


3. EQUILIBRIUM FLOW THEORY FOR AN IDEAL DISSOCIATING GAS 


3.1. Requirements for equilibrium flow 


We consider now the flow of the ideal dissociating gas of §2 under 
conditions in which to a good approximation an equilibrium state is 
maintained at all points. ‘This means (i) that the state of each portion of 
fluid changes slowly compared with the relaxation frequencies of all the 
processes involved in maintaining equilibrium, and (ii) that there is no 
transport of momentum or heat, or significant interdiffusion of atoms and 
molecules, across surfaces moving with the mass velocity of the fluid. 
(In Part II we shall consider flows in which restriction (i1) is removed but (i) 
is retained, and in Part III flows in which both are removed.) 

It is well known that the equations governing such reversibly-adiabatic 
gas motions do not always possess continuous solutions, but that solutions 
involving where necessary a discontinuity or ‘shock wave’, across which 
irreversible changes of state occur but mass, momentum and energy are 
conserved, can always be set up, and these correspond well with what is 
observed. Energy dissipation occurs within shock waves, but the time spent 
by any fluid particles within the region of dissipation is not great compared 
with the relaxation times just mentioned (except for very weak shock waves 
indeed), and, whenever the distance travelled during such a time is small 
compared with the basic length scales of the flow, the solutions which 
treat the shock wave as a discontinuity are useful.. These solutions are 
considered in this section as well as the continuous ones. Indeed, since the 
first change experienced by the air as a supersonic projectile approaches 
is passage through a shock wave, the conditions governing this change will 
be studied first. 
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3.2. Conditions behind a shock wave 

In a frame of reference in which the shock wave is stationary, the 
Rankine Hugoniot equations of conservation of mass, momentum and 
energy across it are 


Po = (27) 

= = (28) 

Pot Pot = Unb (29) 
ig + = + (30) 


where suffix 0 denotes conditions ahead of the shock, suthx 1 conditions 
behind it, and 7, and v, denote velocity components normal and tangential 
to it. 

With atmospheric temperatures ahead of the shock, dissociation can be 
appreciable only if 7, is large —so large, in fact, that the terms py in (29) 
and 7, in (30) can be neglected. ‘Vhus, figures 1 and 2 show that Cancion 
is not very important unless 7, > 0-25, so that by (30) i, +472) > 0:25. 
But 7, < 0-02 if the temperature 7, ahead of the shock is less than 0-006 
(this fraction of 7, is 350 K tor oxygen and 640 K for nitrogen), so that 
in this case % is already small compared with $72, for the rather small 
amount of dissociation implied by the selected value of 7,, and becomes 
quite negligible for larger amounts. ‘The ratio of py to pp Vp is still smaller. 

Accordingly, the ‘strong-shock approximation’ of neglecting 7, and pp 
will here be made, and (27), (29) and (30) written as 


” 


Po no = Pi = Pit Uno = (31) 
It follows from (31) by elimination of v,,, and v,, sia the density-ratio at 
the shock wave is 
—-=,— -l= +1 = K, (32) 
where A is one plus the ratio of the total internal energy of the gas to the 
energy of the translational motion in any one particular direction. Figure 3 
gives graphs of 


as a function of « for different values of log, (1/p,); since interpolation 
with respect to log, (1/p,) is linear for fixed « (that is, at any horizontal 
level), only two curves have been drawn. 

In terms of K, the pressure and enthalpy behind the shock are given 


by (31) as 
1 2 . 1 a | a 
i = ga) (34) 


Note that over 98°, of the kinetic energy of the fluid approaching the shock 
wave 1s converted directly into enthalpy by it 
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To determine conditions behind a shock wave given the density and 
velocity ahead of it, one may first use (34) to give 7, correct to 1 or 2",. 
Assuming a rough value of A of say 10 or 12, one then determines 
from figure 1 the value of « corresponding to the enthalpy 7, and to a 
density Kpy, and hence from figure 3 (using this « and the same value for 
the density) an improved approximation to K. If the first guess of K was 
badly out the process may then be repeated to improve the approximation, 
but otherwise this is unnecessary. 

For example, if in the units here used v,9 = 1:2 and py = 2x 10“, then 
i, = 0-72 and, to a rough approximation (taking K = 12), p, = 2:4 10-® 
and log, (1/p,) = 5-62. From figure 1 for these values « = 0-42, whence 
from figure 3 to a better approximation K = 13-2, giving p, = 2°64 x 10-° 
and log,9(1/p;) = 5°58. If this new value is used instead of 5-62, the values 
obtained for x and K are unchanged in the last figure quoted. 


6 J 6 $ 10 i 13 14 IS 16 


12 
K= 2, /Po 
Figure 3. Density-ratio K at a shock wave, as a function of p,; and a, the density 
and degree of dissociation behind it. Linear interpolation with respect to 
logjo(p q/p1) at any horizontal level is accurate. 

Although we shall not use them in this paper, it is valuable in more 
accurate calculations to have available the second approximation to the 
shock wave equations (32) and (34). This takes into account the terms py 
and 7) in (29) and (30) but neglects their squares. If the gas ahead of the 
shock is taken as an ideal diatomic gas with velocity of sound a, = 1/(7p9/5pp), 
then these second approximations are 


Po HK-1) vo)’ Potro K 7K(K-1) 
2(6K—1) | a, 


The effects on p, and 7, are appreciable (around 5°, even for v,,9/a) = 10) 
but that on p, is very small (less than 0-4°%, under the same condition). 


a 
> 
ate 


14 M. 7. Lighthill 


3.3. Effects of non-uniformity of shock wave strength 

Any projectile travelling at speed V must have a shock wave ahead of it 
which is normal to the direction of motion at one point at least (since an 
ideally sharp apex or edge in front of the projectile could not be used in 
practice). At this point v,,, = V and hence 2, is very nearly }V*. It follows 
from figure 1 that about 5°, dissociation will occur behind this point on 
the shock wave when V = 0:7, 25% when V = 1-0,:'50°, when V = 1:3, 
75°, when V = 1-5 and 95°, when V = 1-6, the unit for V being of 
course v,. ‘These are only rough values; the variation: with density of 
the dissociation produced can be seen from the figure. 

At other points on the shock wave, where it is inclined at an angle 7 
to the direction of motion, v,,. = V sin», and the enthalpy behind the shock 
wave is reduced. ‘Thus, there is less energy dissipation in these weaker 
portions of shock wave, where accordingly the entropy increase is less. 
‘The entropy gradient behind the shock wave is easily inferred from (34). 
If dS, is the difference in entropy behind elements of shock wave with a 
difference dv, in upstream normal velocity, then 


T, dS, = di,- 
Pi 


v2 dK l 1 


the terms due to the variability of A cancelling out. 

This variation of entropy behind the shock leads to the presence of 
vorticity. At a point P where the two ‘ principal curvatures’ of the shock 
surface are «,, and x, (with curvatures taken positive when they are convex to 
the oncoming flow), let suffix a signify components parallel to the line of 
curvature through P which has curvature «,, and suffix 6 similarly. Then 
at a point on the shock surface near P, whose separation from P has 
components 4x,,, dx,, the velocities differ from those at P by the amounts 


But the equation of momentum in the continuous flow behind the shock 
wave may be written in terms of the vorticity w = curl v as 

2,1 9 
wav+V(4v?)+ - Vp = 0, (37) 
p 
where the first two terms are a familiar form of the acceleration of a fluid 
element in steady flow. ‘This, combined with the energy equation 


+7 = constant, (38) 


gives wav a Vinx Vp = TVS. (39) 
p 


& 
q 
J 
: 


Dynamics of a dissociating gas 15 


Now by (36) the vorticity component normal to the shock surface vanishes ; 
thus 


=0. (40) 


Hence, resolving (39) in the a and b directions (assumed such that the a, b 
and n directions form a right-handed system), 


0S, aS, 
a 


From (41), (35), (36) and (27), 
1 T dS, Ono 


(42) 


‘Thus, the vorticity component in either of the two principal directions of 
curvature of the shock surface is equal to the curvature in the direction at 
right angles, times the tangential velocity component, also in that direction, 
times a factor (K—1)?/K, which for a dissociating gas is rather large 


‘(figure 3), 


Results like this are usually stated for two-dimensional or axisymmetrical 
tlow; these cases are simpler because the velocity component in one of 
the two principal directions of curvature vanishes. It may be useful, 
however, to have a more general type of result on record. Equations (42) 
and (43) are also completely general with respect to the thermodynamic 
properties of the gas; only the ‘strong-shock’ approximation has been 
used in deriving them. 


3.4. Properties of isentropic changes 

In the type of flow described in § 3.1, fluid particles undergo isentropic 
changes as they move along streamlines behind the front shock wave (at 
least, until they encounter another shock wave), although the entropy S 
has different values on different streamlines. ‘The character of the flow 
changes when the fluid speed exceeds the local speed of sound, whose 
square is equal to the derivative of pressure with.respect to density at 
constant entropy, even when propagation occurs in a region of varying 
entropy. Both these facts make it necessary to study isentropic changes 
in detail. 

The general character of these changes can be seen from figure 2, two 
points being particularly worthy of notice. The change in? for a given 
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change in pressure is smaller than in most gases. In fact, at constant 
entropy 

dlogp idp ip K+1° 
In the region of dissociation this is of the order of 3, only half of its value 
for an ideal diatomic gas. Hence changes in enthalpy, and so also of the 
fluid speed 


dlogi p 2 (44) 


v= 1/(V?-22), (45) 
with changes in pressure, are reduced when dissociation or reassociation 
is occurring. ‘This is because the translational energy of the molecules is 
a reduced fraction of the total enthalpy. 

Again, the change in logp with logp at constant entropy is slower for 
intermediate x than for either an ideal diatomic gas or a gas with constant 
adiabatic index $ (the limit of our ideal gas as x —- 0), or a monatomic gas 
(x-+ 1). In fact the adiabatic index 

(46) 
dlogp pip 
from which the speed of sound a can be inferred, falls to around 1-2 in the 
region of dissociation. For by putting dS = 0 in (25), and using (23) with 
p expressed in terms of 7 and x. we can show that 


2+3a?—a3 , 


feat 


7 
(47) 
a(1—«) 


From table +, which gives y—1 for different values of x, it is clear that, 
for the values of 7 of interest here (say, between 0-04 and 0-10), (y—1) 
falls well below both its extreme values for all the intermediate values 
of x which are quoted. 


| % 0 0-1 ()-2 0-3 0-4 0-5 
y—1 1 2T-—20:5T? 2T—11T* 2T—-8:-2T? 27 +-7:2T* 2T--7-0T? 
3 1+33-7T* 1+ 24-3T? 1 -—20T? 1--187? | 
y—1 2T-7:°5T° 2T—8-8T? 2T -11-8T? 2T 2 
1+17-5T* 1 ---18-6T? 1--22-5T? 3 
Table 4+ 


This behaviour of y is not unrelated to the behaviour of K shown in 

figure 3. In fact, it follows from (44) that in a change at constant entropy, 
dlogp K+1dlogp—dlogi K+1 dlogp—dlog(}K+}) 

“dlogp K-1 dlogp K-1 dlogp | 

_ K+1—dKk dlogp (48) 


K-1 
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which may be compared with the relationship y = (AK+1)/(A—1) for a 
gas with constant specific heats. ‘The actual expression for the 
isentropic derivative dK/dlogp is algebraically complicated but it can 
be seen from figures 2 and 3 to be numerically small for « < 0-9, so that 
the fall of y to a minimum for intermediate « proceeds in parallel with the 
rise of K to a maximum for these «. 
The square of the local Mach number immediately behind a shock wave, 
by (46) and (34), is 
vite,  _ 1+ K*coty 
since @,9 = Vsiny and v= V cosy. For a dissociating gas, then, this 
number is small (around 0-07) behind the normal part of the shock wave 
(where 7 = $7) but rises to unity where 
(y(K 
7= cor (50) 
which for y = 1:2 and K = 12 (say) is 73-8. Although the sonic line 
starts on a portion of shock wave inclined rather steeply to the flow, where 
the pressure is 


(+9) 


K 
it ends at a point farther round the body, where the pressure (estimated 
from (44) by taking A constant) is about 


1\{ 1 y 


For y = 1-2and K = 12, the pressure (51) is 0-845p,) V? and the pressure (52) 
is 0-540p, V2. 

| 

3.5. Flow about bluff bodies: approximations of Newtonian type — ' 

‘The aerodynamic problems we shall consider will be concerned mostly 
with flow about bluff bodies (with the sphere as the typical body), partly 
because it seems clear that the flow about a projectile must at least start 
by negotiating a more or less bluff nose. 

The theories of flows about bluff bodies at very high Mach number 
which are mentioned in the literature have been related by most authors 
in a somewhat far-fetched manner to a discussion by Newton (1687). 
Having no empirically-based knowledge of gaseous structure, Newton 
postulated a gas which we should now describe as having negligible thermal 
motions and very large mean free path and whose molecules collide 
inelastically with a solid surface. This gives a surface pressure of 
py V?sin?y, where xy is the angle between the surface and the direction 
of motion; all the momentum flux normal to a surface element is communi- 
cated directly to the surface. The reflection condition is satisfied if particles 
hitting the surface are adsorbed and later evaporate at the temperature of 
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the surface with energy small compared with V?sin*y. Cooled bodies may 
therefore satisfy the Newtonian assumptions at Mach numbers so large 
that thermal motions in the undisturbed atmosphere are negligible, but 
only if the density is so small that the mean free path is comparable with 
the body size. ‘The last restriction is a severe one, and applies only at 
extreme altitudes, where the drag and heat transfer are not very important. 

In the more interesting class of flows in which the mean free path is 
small compared with the size of the body, no plausible reason why the 
Newtonian raretied-gas formula might be applicable has ever been advanced, 
but the prestige of Newton’s name has been regarded by some scientists 
as an adequate substitute for a theory. (Newton actually gave a different 
argument, and a different answer, for denser fluids; although these were 
incompatible with our present knowledge of the mechanics of fluids, it is 
significant that he did not himself expect the rarefied-gas formula to be 
applicable. ) 

At the front stagnation point the pressure is equal to the stagnation 
pressure appropriate to the conditions behind the normal part of the shock 
‘To a good approximation, by (49), the flow here may be regarded as 
incompressible, and so, by (34) with v,. = V, this stagnation pressure is 


1 


a value only slightly below the ‘Newtonian’ value py) V*. Lees (1955) 
suggests an empirical law 
p= (1 V*sin*y (54) 
for the surface pressure distribution, in which the Newtonian formula is 
multiplied by a constant factor to make it accurate at the stagnation point, 
and in this way obtains good agreement with experimental pressure distri- 
butions on spheres, and cones of 40° semi-angle with spherical tips, in air 
flows at M/ = 5-8, as observed by Oliver (1956). Such a formula could 
reasonably be used for extrapolation to when dissociation is present only 
if it had some plausible theoretical basis. This seems to be wanting. One 
can hardly regard the flow normal to a surface element (with velocity 
I” sin y) as independent of the flow tangential to it, which completely alters 
the streamline pattern on which the calculation of pressure recovery in (53) 
is based. as well as introducing important centrifugal pressure gradients. 
A quite different approach was set out in a paper by Ivey, Klunker & 
Bowen (1948) and anticipated in a brief discussion by Busemann (1933). 
The approach has been used and further discussed by Grimminger, 
Williams & Young (1950), who call it ‘the hypersonic approximation’, 
and by Van Dyke (1954), who calls it the ‘ Newtonian-plus-centrifugal ’ 
theory. More recently, Chester (1956) and Freeman (1956) have both 
reconsidered the approach and pushed it to a second approximation. 
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‘his approach first assumes the Mach number so great that the strong- 
shock approximation can be used (as we did already in $3.2) and then 
approximates further by taking AK large. ‘This may be regarded as beginning 
an expansion in inverse powers of A, or (speaking more accurately) a 
successive-approximation procedure in which the first approximation is 
the solution with A= «x. Since A, as we have seen, is large for a 
dissociating gas, such a procedure is particularly attractive when dissociation 
occurs. But, as we shall see, it leads to difficulties for certain body shapes, 
principally due to uncertainty about the nature of the limiting flow for 
K = 

Since for K = « the density behind the shock wave is an infinite multiple 
of the density ahead of it, while the velocity component tangential to the 
shock wave is unchanged, the streamtube area might be expected to become 
zero on passage through the oblique parts of the shock wave. In this case 
the shock wave would be wrapped round the body, and all the flow between 
the shock wave and the body would be compressed into a layer of negligible 
thickness. (Note that here viscosity and allied effects are neglected; their 
importance can be gauged better when a second approximation has 
vielded an estimate of the thickness of this ‘compressed-air cap ’.) 

Now, since by (34) for K = oo the pressure immediately behind a shock 
wave is py) V?sin?y, while, if the shock wave is wrapped closely round the 
body, 7 = x at any point, it might be thought that the ‘ Newtonian’ pressure 
distribution py V*sin?y would be recovered in the limit as K-> «©. This 
neglects, however, the centrifugal pressure change across the layer, in which 
a large mass-flow moves with large velocity along streamlines with curvature 
equal to that of the body. ‘The true surface pressure for K = © is 
accordingly less than the Newtonian value. ‘This calculation is briefly 
illustrated for a body of revolution in axisymmetric flow in §3.6 (see 
Grimminger, Williams & Young (1950) for a more general discussion) ; 
some further details found by Freeman (1956) for this case are also described. 


6. Higher approximations: summary of Freeman’s analysis for bodies of 
revolution 
If y signifies distance from the axis of symmetry, and (x,y) are the 
usual boundary-layer coordinates (x = distance along meridian section of 
surface, y = distance normal to surface), and the relationship between r 
and x on the surface is r = R(x), and y is a Stokes stream function, whose 
value at each point is (27)~! times the mass flow across a disc (coaxial with 
the body) with the point on its rim, then in the uniform flow upstream of 
the shock wave ¢% = }p)Vr?, and so on the shock wave (which coincides 
with the body surface for K = 2) % = }p)V R*(x). ‘The angle x(x) between 
the tangent plane and the axis is sin“'R’(x). But the equation of momentum 
normal to the surface gives 


d ap é 


B2 


L 
t 
: 
t 
n 
iS 
3) 
= 
») 
+) 
1S 
Vig 
ly 
ne 
tv 
TS 
9 
) 
ts. 
3, 
3). 
er, 
’ 
al 


20 M. F. Lighthill 


since —dy/dx is the curvature of the streamlines and since ¢cys/éy = pir. 
Here, inertia of flow normal to the surface is neglected, because the flow 
becomes purely tangential as K  «. But the pressure at the shock wave 
was shown to be py, V*sin?y(x); hence, the pressure at an arbitrary point 
behind it (that is, for any value of 4 < 3p) VR(x)) is 


dx/dx 
R(x) J, 


The surface pressure is obtained from (56) by putting 4% = 0. 

One still has the problem of selecting the value of v as a function of u& 
for given x to be put in the integral. Here the properties of isentropic flow 
along streamlines are used. ‘The streamline specified by % passed through 
the shock wave at the point x = x,, say, where % = }p)VR*(x,). At this 
point, by (34) with K = x, 71 = }V*sin*y(x,). Now, equation (44) shows 
that as K -» « the change in 7 with changes in p becomes negligibly slow, 
so that in this limit we may take 

v = (V?—22) = Vos x(x) (57) 


all along the streamline in question. ‘Then (56), with its variable of 
integration changed from : to x,, becomes 


p = po V*sin?y + v dp. (56) 


whence the surface pressure in the form given by Ivey, Klunker & Bowen 
(1948) is derived by putting x, = 0. 

‘The approximations involved in this formula are as follows. 

(i) ‘The shock wave angle 7» has been taken equal to the body surface 
angle x, although for convex bodies 7 > x. 

(ii) The streamline curvature has been taken equal to that of the 
body, although it will be less for convex bodies. 

(iii) The increase in velocity along a streamline as the pressure falls 
has been neglected. 

(iv) ‘The terms in 1/K and 1/K? in the shock equations (34) have 
been neglected. 

(v) The pressure gradient normal to the wall has been taken equal 
to the (centrifugal) pressure gradient normal to the streamlines, 
although since the streamlines are not exactly parallel to the wall 
it must contain a component of the streamwise pressure gradient 
as well. 

Of these, approximations (i) and (ii) tend to underestimate the surface 
pressure, and (iii) and (iv) to overestimate it. Approximation (v) can work 
either way, but is probably the least important of the five. Since the 
experiments of Oliver (1956), although performed at Mach numbers at 
which the strong-shock approximation is in error by 15°,, and under 
conditions in which air behaves like an ideal diatomic gas, with K = 6 
(obviously too low a value for the approximations to be good), give surface 
pressures consistently greater than equation (58) predicts, and in fact 
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close to Lees’s equation (54), we may infer that approximations (i) and (ii) 
are more seriously in error than (iii) and (iv). Probably the order of 
importance of the errors involved in the different approximations is that 
given in the list. 

Attempts have sometimes been made (e.g. by Grimminger et a/. 1950) 
to produce better agreement with experiment by choosing a different 
relation between wv and y% to insert in (56), but this clearly requires v to 
decrease along streamlines, although the pressure is falling. Viscous action 
could explain this, but hardly at the Reynolds numbers of the experiments. 
Probably approximations (i) and (i1) are a much more important source of 
error at fairly low values of A. All the sources of error will, however, 
decrease as K becomes greater. 

In addition to the experimental disagreement, the theory contains 
within itself signs of its own limitations. For certain shapes, the surface 
pressure falls to zero for some x. ‘Thus, for a sphere of radius a, 
R(x) = asin(x/a) and (x) = }+—.x a and the surface pressure given by 
p = py — 4sin*(x/a)}, (59) 
which vanishes at v/a = 47, that is, at 60 from the front stagnation point. 
Ivey, Klunker & Bowen (1948) suggest that the flow separates from the 
surface at this point, leaving a dead-air region near the surface in which the 
pressure coefficient is negligible. They calculate the drag of a sphere on 
this basis as two-thirds of the ‘ Newtonian’ value. 

The present author does not believe that the flow separates from the 
surface at this point, even for very large A, but rather that the shock wave 
separates from the surface there. ‘The assumption of very small streamtube 
area then breaks down, because in hypersonic flow enormous increases in 
streamtube area are possible when the pressure has fallen by a large factor 
(and it must ultimately fall to at least its upstream value). ‘The whole 
approximation thus becomes non-uniformly valid as the point of shock 
wave separation is approached. 

This hypothesis is borne out by the behaviour of the second approxi- 
mation (Freeman 1956) to the Ivey, Klunker & Bowen theory. Freeman 
first finds y as a function of x along each streamline % = constant by 


integrating po ("dys 
60 
in which he takes @ as given by (57), and, since y is so near 1 for K large 
($ 3.4), 
+ pf = (61) 


V*sin®y(x,)’ 
with p given by (58). Here, K(x,) stands for the value of K at a shock whose 
angle of inclination to the oncoming flow is y(x,). Hence, the equation of 
the streamline is 


R(x,) R’(x, dx, 
K(x,)cos sin®x(x) + | R'(x,)cos x(x) ds} 
(62) 
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where x, as before is the position where the streamline passes through the 
shock wave. Note that the distance of the streamline from the surface is 
of order K~'; however, where the pressure p begins to be small, the 
coefficient of A~! starts to increase rapidly. 

The integral (62) converges at its lower limit because R(x,) tends to 0 
as well as cos y(x,) as x, ~ 0. Otherwise (as happens in the corresponding 
theory (Chester 1956, Freeman 1956) for two-dimensional flow) one would 
have to use a better approximation to v than (57) for small x, (i.e. very 
near the body). 

Freeman next obtains the equation y = v(x) of the shock wave, by 
simply replacing the upper limit in (62) by x, so that it becomes the value 
of y for the streamline which passes through the shock wave at the point x 
itself. This expression again is of order K~! but becomes large as the 
point of shock separation is approached. As x — 0, on the other hand, it 
is easy to show that 

a/ K, (63) 


where K = A(0) stands for its value at the normal part of the shock wave 
and a is the radius of curvature of the body at the nose. (It is interesting 
that this stand-off distance of the shock wave ahead of the body is exactly 
twice the value obtained by the very crude approximation of assuming 
irrotational flow behind the shock wave. Equation (63) is closer to the 
experimental results. A more accurate value than either will be found 
in §3.9). Strictly, the assumptions of the theory break down near the 
stagnation point, since the flow velocities are not even approximately 
tangential there, but Freeman (1956) is able to take approximate account 
of the normal components of velocity in this region, and to show that they 
do not affect the equation y = y,(x) of the shock wave to the approximation 
used. This analysis shows also that in this region, for K large, the tangential 
velocity component v, falls linearly with y from its value I”, a on the shock 
wave to zero on the body, a striking result which would restrict greatly the 
importance of viscosity if it were true, and which also will be critically 
evaluated in $3.9. 

He then goes on to obtain the surface pressure distribution to a second 
approximation, in which all the effects numbered (i) to (v) above are taken 
into account, by use of the shock shape and streamline shapes which have 
just been obtained, and an increase of velocity along a streamline which 
corresponds to the pressure fall as calculated on the first approximation. 
To the second approximation thus obtained, the pressure starts at the 
correct stagnation-point value (53), which is below the first approximation, 
and for a sphere (for which numerical details are obtained) begins to rise 
above the first approximation for 6 > 20°, and has already diverged very 
greatly from its first approximation when 6 = 45 . 

It is clear from this that the whole process diverges near the point of 
shock separation, largely because the assumption of small streamtube area 
behind the shock wave starts to break down when the pressure has fallen 


‘ 


Dynamics of a dissociating gas 23 


by a factor at all comparable with the shock-density ratio K. (‘The breakdown 
occurs even sooner in two-dimensional flow.) 

Chester (1956) gets similar results to Freeman in a method of successive 
approximation in which the shape of the shock wave is taken as given and 
the shape of the body progressively determined. He uses a perfect gas 
with constant specific heats but pushes the calculations to one degree of 
approximation higher than those of Freeman. 


3.7. Bodies of revolution with shock wave separation at the base 


It must now be noticed that many bodies of revolution terminate in a 
flat base in such a way that no shock wave separation occurs until the base 
is reached—either because the angle y between the surface and the stream 
remains fairly large ahead of the base, or because as y falls the curvature 
(—dy/dx), and with it the centrifugal pressure drop, falls also. 

For example, on a paraboloid of revolution whose equation in cylindrical 
polar coordinates (r,z) is r? = 2az (here, as in $3.6, a is the radius of 
curvature at the nose), the pressure given by the Ivey, Klunker & Bowen 
formula nowhere falls to zero. For, in cylindrical coordinates, that formula 
(§ 3.6) may be written 

1 d*z/dr*  (dz/dr) dr 
po V2 {1+ (dz/dr)?}?" 

which for the paraboloid is easily evaluated as 

po V2 2\14+ (r/a){1+(r/a)?}?? 
Note that this is always more than half the simple ‘ Newtonian’ value. At 
r = a, where y = 45° (and the meridian radius of curvature has increased 
o 2-83a), p/p) V? has fallen to 0-406. ‘The approximation is of doubtful 
value at pressures much below this, but there is not the same catastrophic 
failure of the theory as there was near the position on the sphere where the 
first approximation to p vanished. 

Another example easily worked out by equation (64) is the body r? = 3622, 
for which the radius curvature falls from infinity at the nose to a minimum 
of 1-306 at r = 0-707b before rising again as r increases further (for example, 
it is 1-415 when r = 6 and y = 45°). The surface pressure distribution is 

1 

giving P/Po V2 = 0-354 ies r=b. Note that the pressure where y = 45° 
is lower in this case than for the paraboloid, whose radius of curvature 
relative to the local value of r is twice as great at that point. The pressure 
at x = 45° on a sphere is lower than either (p/p) V? = 0-333), although the 
ratio of the local radius of curvature to r has the same value v2 as for the 
cubic shape. This is because the flow velocities at y = 45° are smaller for 
the cubic shape, for which less of the fluid has passed through the more/ 
oblique parts of the shock. 


(64) 
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For most other shapes, numerical integration is necessary in (64). One 
can get useful analytical expressions, however, for the purpose of seeing 
quickly how various kinds of variation of shape affect the pressure distri- 
bution, by studying bodies in which cos y is a simple function of r, using 
the form 


(67) 


of the Ivey, Klunker & Bowen formula. ‘The general shape of the body 
is fairly obvious from the form of cos y as a function of r, and can be computed 
quickly from the formula dz/dr = cot x. 

One may use (67), for example, to discuss different ways of rounding 
the nose of a cone of given semi-angle yy. On the cone itself p/p, V? will 
be sin®y,. If a spherical nose is used the pressure just before the conical 
part will be 1—4cos*y, to the first approximation. ‘The sudden rise in 
pressure as the conical part is reached, although it will doubtless be smoothed 
out and reduced in magnitude a bit in the real flow (because the centrifugal 
effect does not cease instantly when the fluid reaches the conical part, and 
because the effect is overestimated in the Ivey, Klunker & Bowen formula), 
may cause undesirable boundary layer behaviour (instability, etc.). With 
a paraboloidal nose the pressure falls instead to 


I(sin® yy + sin* yg sec sinh™'cot XQ), 


for example (if yy) = 45°) to 0-406 instead of 0-333, but the sudden rise to 0-5 
may again be harmful. Instead, one can let the radius of curvature become 
infinite where the nose joins on the cone. Choosing 


COS X = SEC Xo» (63) 


which gives y = yy and dy/dr = 0 when r = 2acos yo, we tind from (67) 


that 
4fr\?  35/r\3 3/r\t 
(;) SEC Xo— } SEC*Xo- (69) 


This is continuous with p/p, V? = sin?y, at r = 2acos x», but falls to a 
minimum of | — 1-055 cos*y, at r = 1:735a cos xo, for example to a minimum 
of 0-473 if y, = 45°. Thus some slight pressure rise still occurs but it is 
much smaller and more gradual. The need for some pressure rise 
immediately before an accurately conical portion is easily seen by 
differentiating (67). 

Conversely, equation (67) may be inverted, to give the shape of body 
with a given first approximation to the pressure distribution, as 


| (1 dr, 


(70) 


cosy = 
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(which must be used with dz/dr = cot x if the body shape is to be plotted). 
Here p, is the pressure (to the first approximation) at the point r = r, on 
the surface. 

For example, if one takes the pressure to be its value on a sphere of 
radius a up to a certain value of r, and thereafter constant and equal to that 
on a cone of semi-angle yp, in other words if 


‘ 
pi 4 
r < ahv3 cos xo), 
Po V 2 3 a ( 2 Xo) (71) 
= sin* yy (r > cos x), | 
then (70) gives 
cosx =ra (r < cos | 
| 
3 \—1/2 (72) 
COS Xo 
= COS xy, 1 + (r > ahv3 cos xp). 
Xo} — 3a? 2 Xo) j 


Sphere 


Figure 4. Transition curves from spherical nose to cones of semi-angles 45° and 30”. 
Bodies of revolution with these meridian sections would have constant surface 
pressure behind the spherical portion, according to the Ivey, Klunker & Bowen 
formula (67). 
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It is to be noted how rapidly . > x» as r increases; the transition from 
sphere to cone is made fast, although without the pressure dropping below 
its final value. ‘The shape is plotted for yy = 45° and y, = 30 in figure 4. 
These are probably rather suitable ‘transition curves’ from the spherical 
to the conical shape, since the fact that the pressure formula overestimates 
the centrifugal effect can only mean that the tendency to pressure rise, 
already eliminated by the shape chosen, is less marked even than was 
assumed in the calculation. 


3.8. Shape of shock wave beyond where it separates from the body 

Both for bodies from which the shock wave separates at the base, and 
tor those from which it is thrown off (as it were) by centrifugal force at a 
point ahead of the base, it is of interest to investigate the shape of the shock 
wave beyond where it separates. “lhe following rather crude approximate 
treatment of this problem is based, as in $$ 3.6 and 3.7, on assuming K large. 

For K = « the shock wave coincides with the body up to the position 
of shock separation, and all the streamlines coincide with both shock and 
body. After shock separation, however, it seems likely that most of the 
streamlines will remain close to the shock rather than to the body, for when 
the shock is oblique the streamtubes passing through it are greatly reduced 
in area and lie closely parallel to it, and, also, all those streamtubes on which 
the pressure has not fallen off greatly from its value at the shock will still 
have a smal! streamtube area. If the shock has not become too oblique to 
the flow, then most of the mass flow inside it must flow along these stream- 
tubes. ‘The streamtubes nearest the body or wake, however, must have a 
low pressure, either because the wake is at a low ‘ base’ pressure or because 
a large streamtube expansion in this region is necessary to fill up the growing 
space between the shock wave and the body. This suggests that the shock 
wave has such a shape that the centrifugal pressure drop, across streamlines 
most of which nearly coincide with it, reduces the pressure from its value 
at the shock to practically zero on the other side of those streamlines. ‘This 


ives 
0 = py V? sin*y + 


dis 


r(x) ~0 
where x 1s now a coordinate measured along the shock instead of along the 
body surface, 7,(x) is the value of 7 at the shock, and the streamline curvature 
has now been taken equal to the shock wave curvature (—dn/dx). Because 
r(x) = sin n, this equation may be written 


sing! + | r,cosy dr,= 0, (74) 


an equation which is ov solved after differentiation with respect to 7 as 


const. = —siNn { 7.0087 dr., (75) 


sin? 
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where suffix zero signifies the position of shock wave separation. (To a 
first approximation, one may take y = y and r, = R(x) on the right-hand 
side of (75), since the shock coincides with the body up to separation.) 
Next, integrating (75) again, we obtain 
— 7?) = (cot 7 — cot H9)sin | r.cos dr, (76) 

as the relation between 7 and r, along the shock. (One may note that this 
is also derivable, less simply, from (70) by putting p, = 0 for r, > ry and 
interpreting y as 7.) 

The shock shape in cylindrical polar coordinates (r, z) is easily deduced 
from (76), because dz,/dr, = coty. Hence 


6" § 


= yo + (77) 


SiN Ho | r,cos dr, 


showing that the approximate shock shape calculated in this way is a cubic. 

The errors in this formula are similar to those enumerated after equation 
(58), except that (i) is absent; no error has been made in the shock wave 
angle in equation (73). However, the error in the streamline curvature 
may be more severe in this case. ‘The streamlines near the body or wake 
are probably considerably more curved than the shock wave, indicating 
that (77) probably overestimates the curvature needed in the shock wave. 
One may hope, however, that their contribution to {vd is only moderate 
because neither v nor the change in & across those streamlines can be very 
large. At present it does not seem clear how the theory beyond shock 
separation can be extended to a second approximation. However, in cases 
where separation occurs at the base, one may use Freeman’s calculation 
ot the second approximation to the shock shape up to separation, and take 
39s To No IN (77) as the coordinates and angle of the shock wave at separation, 
and thus considerably improve the accuracy of (77) beyond separation. 

When separation occurs before the base, the theory is less satisfactory 
because for any finite K a substantial departure of the shock wave from 
the body surface occurs before the theoretical position of separation, and 
Freeman’s second approximation is of no value for calculating the later 
stages of this process. In these cases, the first approximation to the shock 
shape is put forward here principally because it might serve as the foundation 
of a future method of successive approximation, based (as is necessary for 
success) on a uniformly valid first approximation. 

Figure 5 shows the shock wave shapes for separation from a spherical 
cap subtending a half-angle at the centre of 10°, 20, 30°, 40°, 50° and 60°; 
for convenience of comparison they are all superimposed on a single figure, 
and the values for K = « only have been plotted. ‘The first three are 
likely to have more value than the others; for finite K these may well be 
close to the true shock shape if the body is taken to be a portion of a sphere 
concentric with the spherical part of the shock but of a smaller radius (to 
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be determined in §3.9; the position of such an inner sphere for K = 15, 
12, 9, 6 is noted on the figure). The last of the shock shapes shown is a 
very crude first approximation, for large K, to the shape of the shock wave 
on a complete sphere (which separates from the surface at y = 30°). 


Position of 
spherica 
cap when it 

K 


Figure 5. Shapes of shock waves trailing behind spherical caps subtending semi- 
angles 10 , 20 , 30°, 40 , 50°, 60° at their (common) centre. 


3.9. Constant-density approximation to flow near front stagnation point 


An approximation for use in the neighbourhood of the front stagnation 
point (a region important for heat transfer since the fluid temperature is 
greatest there and the boundary layer thinnest) is now obtained without 
any procedure of successively discarding inverse powers of K. It can be 
used to check some aspects of theories based on those procedures. 

‘The present approximation is based instead on the theory of incom- 
pressible flow; thus it neglects variations of density, both (i) along the back 
of the shock wave——this is permissible wherever 7 does not vary much 
from $7, in which case v,,, = V'siny cannot vary enough to produce any 
significant variation of A = p,/p)-—and (ii) along streamlines behind the 
shock wave—this confines us to regions in which the velocity does not 
approach close to the local speed of sound. Both assumptions should be 
adequate in the region behind that part of the shock wave in which 7 varies 
between 80° and 90. 

‘The body shape in the region of interest is taken to be part of a sphere 
of radius a, since bodies with at least a spherical portion of surface near 
the nose are of great interest. For other bodies of revolution, however, 
the theory could reasonably be applied with a as a mean radius of curvature 
of the surface in the region just described. 

‘Taking incompressible flow behind the shock wave, and A constant in 
the boundary conditions at the shock wave, we shall show that all the 
equations of motion and boundary conditions can be satistied by assuming 
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a spherical shock of radius c >a. We use spherical polar coordinates 
(r, 8, A) with the line @ = 0 pointing upstream, and a Stokes stream function & 
such that a 
= (rsin @)v,, 
where the suffixes on the v’s signify components. 
On the spherical shock wave, by (27) and (28), 


—(r* sin @)v,, (78) 


v, = Vsin8, (79) 
and these conditions are satisfied if 
Vc? sin?6 Cys 
= = Vesin? =c 80 
yb Vesin?@ onr=c. (80) 


In axisymmetrical flow the only component of vorticity is 
OF rsin@er? 06 
But on the shock wave, by (43), 
(K—1)? V sin 
K 
Further, in axisymmetrical flow, the ratio w,/(rsin@) remains constant 
along each streamline, since the intensity of each vortex ring changes 
during convection in proportion to its radius. But since 
—(K-1PV 
rsin6 Ke? 
on the shock wave, and all the streamlines pass through the shock wave, 
the ratio must have this constant value everywhere. 
Hence, by (81), 


-_1)2 Vr2sin2 
ep | sind é (cosee (K — 1)? Vr? sin?6 (84) 


(83) 


everywhere. ‘This equation will now be solved subject to the boundary 
conditions (80), and it will be shown that the dividing streamline is a sphere 
of radius r=a<c. (Note that no boundary condition on the pressure 
need be applied at the shock wave; the equation of motion determines the 
pressure once the velocity field is known, and the condition (43) on the 
vorticity was obtained as the condition that the boundary condition on the 
pressure is compatible with the equation of motion.) 

The solution of (84) under the boundary conditions (80) takes the form 


= (85) 
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The function ‘(r) satisfying these conditions is easily found to be 

Tir) = 0K’ 3(K 1) (:) —5K(K—4) @) 


| 
(88) 


The solution given by (85) and (88) shows that the streamline % = 0 
divides at r = a, the greatest root of ‘’(r) = 0 which is less than c. (It can 
be shown that roots of ‘f(r) = 0 with 0 <r <c exist for all K>1.) The 
streamline then continues along the sphere r = a, and hence the flow may 
be regarded as the flow about a solid sphere of radius a. 


Figure 6. Ratio of stand-off distance of shock wave ahead of sphere at radius of sphere, 
as a function of K. 


Particular interest attaches to (c—a)/a, the ‘stand-off distance’ of the 
shock wave as a fraction of the radius of the sphere. This is a function 
of K alone, a graph of which for K > 4 (an inequality which by the physical 
significance of K must be satisfied for all gases) is given in figure 6. 
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One is interested also in the distribution of v, between the shock wave 
and the sphere. By (78), Y’(r) (89) 
V sind Vr 
This is plotted against r/c for a number of values of A in figure 7. The 
flow is seen to be in nearly (though not exactly) uniform shear, as predicted 
by Freeman’s first approximation (§ 3.6), but the velocity at the wall is not 
zero as that approximation would give. Freeman’s second approximation 
to the velocity at the wall (derived from Bernoulli’s equation and the Ivey, 
Klunker & Bowen value of the pressure) is, however, accurate to within 1°, 
for all the A’s shown. Putting this a different way, the pressure on the 
surface derived from (89) is 


= pe § + 0-01)sin?6 > (90) 


that is, the maximum departure of the coefficient of sin?# from the Ivey, 
Klunker & Bowen value § for K > 4 is 0-01. 


0 i j 
0:38 0-96 0:94 0:92 0:90 0.88 0-86 0-84 


Figure 7. Distribution of transverse velocity zy in the nose region of the flow around 
a sphere, plotted for different values of K (the density-ratio across the normal 
part of the shock wave). The coordinate r (the distance from the centre 
of the sphere) takes the value c at the shock wave; the position of the body 
surface (r = a) in each case is shown by hatching. 


‘This agreement in the rate of pressure drop along the surface, near the 
nose, between two very dissimilar theories (the present one, which only 
approximates by taking the density constant, and the Ivey, Klunker & Bowen 
theory, which allows density variation but assumes large A, and hence a 
large reduction in streamtube area behind the shock wave) gives some added 
support to both, and makes one reasonably confident about the general 
view of equilibrium flow around bluff bodies, in the region where pressure 
is not a small fraction of p, V*, presented in this $3. What should be done 
to treat the flow where, as a result of the pressure falling by a large factor, 
the streamtube area becomes large again, remains at present a mystery. 
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SUMMARY 

This paper is a discussion of recent experiments in shock-wave 
refraction which have clarified a special type of shock outflow 
process appearing to have relevance to other shock interactions, 
and notably to shock reflection from an oblique wall. For certain 
incident shock strengths and angles of incidence z, the air/methane 
refraction problem simulates closely the situation in the trouble- 
some range of the reflection problem, in which « lies between the 
value «, at which the theoretical solutions terminate and the 
value %) that marks the onset of Mach reflection, and in which 
the flow deflections cannot be reconciled with theoretically 
permissible reflected shock strengths. In the analogous refraction 
cases, the reflected shock is observed to increase in strength along 
its length to a maximum value at the intersection point, and to be 
followed by a subsonic rarefaction zone which also increases in 
severity near the intersection. In fact, this zone appears to 
coalesce into a subsonic discontinuity, just at the intersection 
point—a feature which would contradict one of the basic assump- 
tions of the regular reflection and refraction theories. Other 
refraction experiments suggest that a similar process is relevant 
to the Mach reflection configuration, and may account for the 
discrepancies in the three-shock theory for weak incident shocks. 


1. INTRODUCTION 
In the branch of fluid mechanics concerned with shock-wave phenomena, 
the problem of the reflection of a shock front from a rigid boundary has 
long been a matter of fundamental interest and concern. ‘This was the 
first shock interaction to be studied experimentally in sufficient detail to 
permit a significant appraisal of the theoretical techniques available for the 
treatment of these non-linear problems. ‘The earliest results from such 
experiments uncovered discrepancies between mathematical prediction and 
observed behaviour that were sufficiently acute to cast doubt on the validity 
of some basic theoretical assumptions. And to this day, despite considerable 
systematic experimental and theoretical work on this subject (Smith 1945; 
Fletcher, Taub & Bleakney 1951; White 1952; von Neumann 1943; 
Polachek & Seeger 1944; Bargmann 1945; Bleakney & ‘Taub 1949), 
certain of these paradoxes have not yet been resolved, and our understanding 
of the reflection problem remains incomplete. 
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The theoretical approach to the problem is discussed in detail in most 
of the references cited above, and has recently been reviewed by Griffith & 
Bleakney (1954). For the purposes of this article it will be sufficient to 
recall the assumptions that underlie the formulation of the problem* (see 
figure 1). 

(a) ‘The interaction of the incident shock with the wall results in a single 
reflected shock travelling away from the wall into the medium behind the 
incident shock. 

(b) Each of the three angular regions of flow formed by this configuration 
of shocks and boundary is uniform; hence the state of the gas changes 
only across the two shocks. 


I 
R 
a 


Figure 1. Regular reflection of a plane shock : incident shock J deflects flow toward 
wall; reflected shock R re-deflects flow parallel to wall. 


(c) Each shock can be treated as in the standard Rankine-Hugoniot 
theory to relate the change in flow velocity to the pressure ratio across it. 

(d) ‘Vhe net deflection of the gas flow by the two shocks is such that 
the flow in the region behind the reflected shock is parallel to the wall. 

(e) ‘The configuration is stationary when expressed in coordinates x/t, 
y/t. 

(f) No energy is lost to the wall during the process. 

Assumption (6), which will require further examination later, was 
introduced to simplify the analysis. Without it, one would need to include 
the appropriate gas-dynamical equations for the two-dimensional flow in 
the angular regions, to be solved simultaneously with the shock equations. 
This assumption has the effect, in the more general case of curved shocks 

* The experiments on which this paper is based were performed in a shock tube 
(cf. Bleakney, Weimer & Fletcher 1949). Unless otherwise noted, all discussions 
and illustrations likewise refer to the shock-tube situation, in which the incident 
shock advances into a gas at rest. In comparing shock tube problems with the 
analogous wind-tunnel or free-flight situations, it should be noted that more than a 
simple coordinate transformation is involved. In the latter cases, boundary layers 
exist on all surfaces exposed to the flow; and it is with these, rather than with the 
surfaces themselves, that the shocks interact. 
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and non-uniform flow fields, of restricting the consideration to a small 
area surrounding the intersection point—small enough that the changes of 
state of the gas in the angular regions are negligible compared with the 
discontinuous changes across the shocks. 

The assumptions listed above define a ‘regular reflection’ process, 
solutions for which can be found only for angles of incidence « smaller 
than some limiting angle «,, which depends on the incident shock strength 
and the gas being considered (see figure 2). For larger values of x, no real 
solutions exist, implying that physically the assumed process cannot occur 
for these angles of incidence. 


Figure 2. Critical angles of incidence in shock reflection: theoretical limit for 
regular reflection «,, experimentally observed onset of Mach reflection a; 
and condition for sonic outflow «,; 7s. inverse pressure ratio across incident 
shock front. 


In most respects, experimental observations support the regular 
reflection theory. When a plane shock wave is caused to impinge on a 
rigid wall at any angle smaller than ~,, a reflected shock does appear, inter- 
secting the incident shock at the wall, of strength and angle of. reflection 
in good agreement with the theoretical predictions (cf. Smith 1945). 

Extension of such experiments to larger angles of incidence reveals 
that beyond some limiting value %, which is very close to (but always 
slightly larger than) the theoretical value x,, the intersection point leaves 
the wall, and occurs instead at some finite distance from it. ‘This inter- 
section is then joined to the wall by a third shock, called the stem, and also 
by a slipstream, or contact surface, which trails behind the pattern, thus 
forming the well-known Mach configuration (see figure 3). The three- 
shock interaction occurring in this pattern has also been studied theoretically. 
This theory presumes that the pressure changes and flow deflections 
experienced by the gas passing through the pattern immediately above the 
‘triple point’ (i.e. the intersection of the incident, reflected, and stem 
shocks) are the same as that for the gas passing immediately below it, although 
the exit velocities may in general be of different magnitude (as evidenced 
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by the slipstream). Again, uniformity of the angular regions of flow is 
assumed, restricting the consideration to the immediate neighbourhood of 
the intersection. 


2. ‘THE SHOCK REFLECTION PARADOX 
Two areas of disagreement with the theories outlined above may be 


demonstrated experimentally. For convenient reference later, they are 
labelled I and II as follows: 


I. The regular reflection process appears to persist for a small range 
of angles of incidence beyond the theoretical limit «,. The inter- 
section point remains at the wall for angles past ~, up to %», at which 
angle the Mach configuration begins. ‘The interval (%)—«,) is 
measurable although small, being about two or three degrees for 
intermediate shock strengths (see figure 2). 

II. Although the theoretical solutions for the Mach configuration 
(i.e. for the strength of the reflected shock at the triple point) are 
tolerable for strong shocks, they become increasingly inadequate 
for weaker incident shock strengths (see figure 4). 


Figure 3. Mach reflection : intersection of reflected shock R with incident shock J 
is joined to wall by stem shock S. Slipstream (contact surface) S\S separates 
two parallel outflows of different velocity, density and entropy. 


Note that both of these problems involve essentially the same difficulty, 
namely the existence of a reflected shock that appears to violate the boundary 
requirements of the flow. 

These discrepancies would not be anomalous if it could be shown that 
one of the assumptions of the theory is physically inadequate in the regions 
of discord, or, alternatively, that the experiments have insufficient resolution 
to reveal the true nature of the interaction in the small region to which the 
localized theory applies. Attempts to confirm such suspicions by shock 
reflection experiments have so far been inconclusive, however, as they are 
hampered by two inherent disadvantages. First, the angular interval of 
interest (%)—,) is so small that, considering the resolution attainable in 
such experiments, it is difficult to make any systematic study of the angular 
dependence of a process within this range. The second, and probably 
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more significant complication is introduced by the extraneous disturbances 
from the physical corner at which the reflection process begins. It is a 
curious feature that the outflow from the shock reflection pattern is subsonic 
for angles of incidence greater than a certain value «,, which is just slightly 
smaller than «, over the entire range of incident shock strength (see figure 2). 
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Figure 4. Shock reflection observations : experimental values agree well with regular 
reflection theory in the range ~a< «,. At a,, the observations leave the 
theoretical curve tangentially, and beyond ap, in the Mach reflection region, 
disagree completely with the three-shock theory. (w, w’ are the angles of 
incidence and reflection appropriate for pseudo-stationary description of the 


problem. For a < %, w = a,w’ = 


‘This means that the state of the region behind the reflected shock, which 
undoubtedly is the most interesting, is determined not only by the reflection 
process, but also by the nature of the boundaries at the ‘corner’ where the 
process first begins. Consequently, although some interesting density 
irregularities are observed in this region both for «, < « < % and for the 
Mach patterns, it is difficult to determine whether these are essential parts 
of the reflection process or merely evidence of the corner interference. 
(In some respects it is tempting to consider that it is the presence of the corner 
which causes troubles I and II, but there is little experimental support for 


this point of view—cf. § 6.) 


3. SHOCK WAVE REFRACTION 

Somewhat less ambiguous information on this subject has recently 
appeared in the course of an experimental study of another type of shock- 
wave interaction, namely the refraction of shock waves at interfaces between 
two gases (see Jahn 1956). This problem, while perhaps one degree more 
complex theoretically than the reflection problem, bears many similarities 
to it; and in some respects it is more tractable experimentally. ‘The theory 
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of shock retraction (Taub 1947; Polachek & Seeger 1951) assumes that a 
pseudo-stationary pattern will be formed, consisting of three shocks or two 
shocks and a Prandtl-Meyer wave (for our present purposes we need 
consider only the former case), and intersecting at a point on the gas 
interface (see figure 5). It is then required that the exit pressures, and the 
flow deflections via the upper and lower paths, be the same. Again, for 
given incident shock strength and gas combination, theoretical solutions 


Figure 5. Regular refraction of a plane shock : shock J, incident on interface between 
two gases O, produces a reflected shock R and a transmitted shock T, the 
outflows from both of which are parallel to the deflected interface D. 


Figure 6. Mach refraction : intersection of reflected shock R and incident shock J is 
joined to interface by stem shock S. Again there is a slipstream SS from 
the intersection point. 


can be obtained only for angles of incidence smaller than some critical 
value «,. Furthermore, it appears experimentally that the assumed con- 
figuration persists, at least locally, beyond «, up to some larger angle x 
at which a ‘ Mach-refraction’ pattern sets in (see figure 6). 

Here too, then, there is a region of disagreement with the theory, akin 
to discrepancy | in the reflection problem. It is not unreasonable to suspect 
that the fundamental difficulty may be the same in both cases. Fortunately, 
the situation is more amenable to experimental study in the refraction case, 
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since the interval (x) — x,) can be much larger for easily realizable conditions. 
Furthermore, in refraction experiments it is possible to change the arrange- 
ment of the physical boundaries to some extent, without necessarily changing 
the nature of the refraction conditions themselves; and thereby the features 
of the pattern which are consequences of the boundary conditions may be 
separated from those which are inherent in the shock interaction. 

The results of refraction experiments pertinent to this problem are 
illustrated in figure 7 (plate 1). ‘These are shock-tube photographs, taken 
through a 5-inch Mach-Zehnder interferometer set in parallel fringe 
adjustment, of the refraction of shock waves with pressure ratio 3-3 at an 
interface between air and methane. Figure 7 (a) (plate 1) shows the con- 
figuration for « <«,. Note that the three shocks are straight, the regions 
between them are uniform, and the intersection point is on the interface. 
Such a pattern, called a‘ regular refraction’, corresponds to the theoretical 
refraction solutions for this problem. Actually, the correspondence is 
encouragingly close; measurements of the shock strengths and angles of 
the configuration, taken from interferograms such as this, have been found 
to be in excellent agreement with the quantitative predictions of the theory. 


Figure 8. Pressure profiles of the reflected wave appearing in figure 7 (4) (plate 1). 
Shock front and following rarefaction zone Z are stronger near the interface 
(dotted profile) than further from it (solid profile). 


However, as the angle of incidence is increased past x, into the critical 
region %, << < 4%, the pattern shown in figure 7(b) is observed. Although 
the shocks still intersect at a point on the interface, this is a more complex 
interaction in which the reflected wave is no longer straight, uniform in 
strength, or ‘flat-topped’. Rather, it is considerably stronger near the 
interface, and less inclined to the flow there than it is further out into the 
field, and it is followed immediately by a very sharp rarefaction region, 
which also increases in severity near the interface. ‘This composite reflected 
signal is thus a form of peaked shock wave (see figure 8). 

The third interferogram in this series (figure 7(c), plate 1) shows the 
pattern for % > %. Here the reflection joins the incident shock wave a 
short distance away from the interface, and a stem shock and slipstream are 
now present in the pattern. Note that again the reflected wave is curved 
and non-uniform, and is followed by a region of ratefaction. 

It is the non-uniform reflected shock and the sharp rarefaction zone 
seen behind it in figure 7(b) (plate 1) and, to a lesser extent, in figure 7 (c) 
(plate 1) which particularly attract our attention in these configurations. 
A similar series of patterns, displaying zones of similar density variation 
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can also be observed in shock-reflection experiments performed at the 
appropriate values of « (see figure 9, plate 2). ‘These zones were noted by 
several earlier reflection investigators, but the unavoidable interference of 
the corner signals largely precluded evaluation of their significance to the 
interaction. ‘lhe contribution of the refraction experiments lies in the 
possibility of eliminating the corner interference in the region of interest, 
so that the observed effects may be attributed directly to the shock interaction. 
It is this technique which permits us to attach significance to the sharp 
rarefaction zones seen in the refraction patterns displayed above, and then, 
because of the near equivalence of the two problems, to regard the corre- 
sponding zones of the reflection configurations in the same light. On this 
basis, all of the remarks in the —— discussion pentem: to both the 
reflection and refraction situations* 


4. THE REFLECTED WAVE 

In addition to the more obvious features of the zone of interest around 
the reflected wave (i.e. the severity and direction of the density gradients 
and the variation of the strength of the reflected front itself along its length), 
two other pieces of information can be extracted from the interferograms. 
‘These are (1) that the region behind the reflected shock is one of subsonic 
flow, and (2) that the strength of the shock at the intersection point is very 
sensitively dependent on the angle of incidence. 

The subsonic character of the flow behind the shock is revealed experi- 
mentally by the catching-up of the sonic signals from the corner (in those 
problems in which they are not exactly compensated to zero strength), or 
it can be deduced directly from measurements of the incident and reflected 
shock strengths and their inclinations to the flow (see Bleakney & Taub 
1949). That this is a subsonic field suggests, of course, that any flow 
deflection processes taking place here will be continuous, relatively wide- 
spread adjustments, rather than the discontinuous, sharply localized 
processes found in supersonic flows. ‘The significance of this point is 
emphasized by the work of Smith (1945), and Bargmann & Montgomery 
(1944), who have shown that a theory which assumes any supersonic 
variations in this region, such as other shocks, or Prandtl-Meyer fans, is 
in even poorer agreement with experiment than the theory postulating 
uniform flow. 

It was mentioned that the strength of the reflected shock front increases 
rapidly along its length to some maximum value at the intersection point. 
This in itself suggests that some unusually severe process is taking place 
at the junction of the two shocks. It is also striking to plot the observed 
strength of the reflected shock at the interface against the angle of incidence, 
from which it appears that for angles immediately beyond «,, the shock 


* The attention of the reader is called to the recent work at Princeton of W. R. 
Smith (1956) in which reflection paradox I has also been examined in the light of a 
related shock interaction—in this case the oblique intersection of two shock waves 
of equal strength. 
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strength is extremely sensitive to «, and rises to a rather surprising value 
just before a (see figures 4 and 10). For angles beyond a, the reflected 
shock strength falls off rapidly with increasing x, suggesting that the 
departure of the intersection point from the wall (or interface) relieves the 
situation which precipitated the sharp increase at x,. ‘That this is indeed 
the case may be seen in figure 7(c) (plate 1). In this well-developed Mach 
pattern, the slipstream from the intersection, which assumes the direction 
of the local flow velocity, is not parallel to the interface, indicating that 
the reflected shock has not had to deflect the flow through so great an angle 
as in the regular configurations. 
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Figure 10. Shock refraction observations for « — a, : strength of reflected shock 
front RS (pressure behind/pressure ahead) vs. angle of incidence «. 


It appears then that we are dealing here with a transition process which 
prevails only in the brief angular interval between the degeneration of the 
regular pattern at ~, and the onset of the Mach configuration at %. This 
transition pattern is somewhat hybrid; it is ‘regular’ in the respect that the 
intersection point occurs on the boundary, but certain features of it, 
apparently associated with the non-uniform reflected shock and the subsonic 
rarefaction field following it, are somehow incompatible with the regular 
theory. 

5. INTERPRETATION OF THE OBSERVED REFLECTED WAVE 

Some insight into the significance of the anomalous reflected wave 
observed above may be gained from examination of other, more familiar 
gas-dynamic situations in which similar flow deflection requirements arise. 
For this purpose it is helpful to recognize that the limit ~,, at which these 
effects first appear, arises as a consequence of the well known ‘maximum 
deflection’ condition of oblique shocks (see, for example, Liepmann & 
Puckett 1952). At «,, the reflected shock is of that optimum strength and 
corresponding inclination to the flow which accomplishes the maximum 
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possible deflection of that flow. Beyond x,, the gas outflow from the 
incident shock is such that no stationary oblique shock can re-deflect it 
parallel to the wall. ‘The situation at the extreme angle in the refraction 
problem, while complicated by the adjacent conditions for the transmitted 
shock, is fundamentally the same. Here there are two shocks (reflected 
and transmitted) standing in two converging uniform streams. Under the 
requirement that the exit pressures be the same, there is just one arrangement 
of strengths and accompanying positions of these shocks by which the 
maximum deflection of the converging flows is accomplished. At ~,, this 
optimum configuration has been established, and beyond it some new process 
must occur. 

An analogous flow deflection limit arises in the more familiar problem 
of the supersonic wedge aerofoil (figure 11, plate 3). Here, -the shock tube 
is used as a transient wind tunnel. ‘The incident shock has passed entirely 
over the model (from right to left), and we are observing the behaviour of the 
aerofoil in the ‘steady’ supersonic flow set up by that shock. (In this case, 
of course, a boundary layer develops along the aerofoil surface, but it is of 
zero thickness at the leading edge and does not fundamentally change the 
nature of the deflection problem.) Over a certain range of conditions, the 
deflection of the flow parallel to the aerofoil surtaces is accomplished abruptly 
by two straight bow shocks, attached to the nose of the body (see figure 11 (a), 
plate 3). (To emphasize the analogy this might be called a ‘regular’ 
deflection.) Beyond a certain limiting body angle (or below a certain flight 
velocity), however, it becomes impossible for any stationary oblique shock 
to deflect the inflow adequately. In this case, the bow shocks are observed 
to become stronger near the leading edge, and less inclined to the flow there. 
These stronger shocks produce subsonic, rather than supersonic, regions 
behind them containing expansions which are most severe behind the 
strongest portions of the shocks (figure 11(4), plate 3). As the limit is 
exceeded further, the bow shocks fuse into one continuous front, detach 
from the leading edge, and advance into the oncoming flow (figure 11 (c), 
plate 3). 

The interpretation of the processes seen in figures 11(4)and 11 (c) (plate 3) 
seems straightforward: whatever part of the necessary flow deflection that 
is not accomplished abruptly through the bow shock is taken care of in a 
continuous fashion in the subsonic field that follows it. Such a subsonic 
field requires a stronger shock ahead of it than the supersonic field of 
figure 11(a) (plate 3), and hence the shock must be more nearly normal to 
the inflow. If the deflection situation is made yet more severe, the region 
behind the shock becomes further subsonic, requiring a still stronger, more 
normal shock. Beyond the condition which produces a bow shock that is 
just normal to the flow at the wedge tip, there is no mechanism for further 
increasing the shock strength, other than for it to advance into the oncoming 
How (thereby increasing its Mach number). ‘The gap thus created between 
the wedge and the bow shock relieves the deflection requirement somewhat, 
and an equilibrium situation is established with the bow shock fixed a finite 
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distance ahead of the wedge (figure 11(c), plate 3). (This equilibrium 
position has been discussed by Laitone (1955).) 

‘The similarities in the nature of the flow deflection requirements, and 
in the observed behaviour of the shocks and the regions behind them, in 
this problem and in the reflection—refraction cases are sufficiently striking 
to permit the inference that the mechanism in all three problems is 
fundamentally the same. In other words, the interpretation applied to 
the three aerofoil patterns in figure 11 (plate 3) should be equally appropriate 
for the regular, transition, and Mach patterns of the reflection-refraction 
problems. Such an analogy implies that the reflection and refraction 
patterns pass through the following sequence of phases as the angle of 
incidence is increased. From normal incidence up to x = «,, the necessary 
re-deflection of the flow parallel to the boundary is accomplished abruptly 
by areflected shock. Beyond ~,, no reflected shock is adequate to accomplish 
this re-deflection entirely, and the remainder of the process is brought about 
continuously in a subsonic rarefaction field that follows the retlected shock. 
Such a subsonic region requires a stronger reflected shock front, which in 
turn must be more normal to the outflow from the incident shock. As « is 
increased further beyond x,, this outflow decreases in velocity, while 
simultaneously the deflection requirement becomes more severe, both of 
which force the reflected shock to be yet more normal to the flow. Conse- 
quently, at some larger angle %, the reflected shock has become just normal 
to the flow, and has no further recourse than to advance into it, much as 
the bow wave detaches from the aerofoil. In so advancing, the reflected 
shock overtakes a portion of the incident shock, thereby forming the 
characteristic Mach pattern. ‘lhe departure of the intersection point 
from the interface relieves the deflection requirement somewhat, so that 
for a given angle of incidence an equilibrium condition is established with 
the intersection point fixed on the incident shock (in the pseudo-stationary 
frame of reference)*. 

Superficially, this conception of the transition reflection process might 
not seem incompatible with the theory outlined earlier; but actually it 
does involve a feature which contradicts one of the basic assumptions. 
The flow deflection adjustment by the subsonic field following the shock 
has been described as a ‘ continuous’ process, as it appears in the interfero- 
grams, to distinguish it from the abrupt change which occurs through the 
shock itself. ‘The only portion of this field which is pertinent to the local 
theory, however, is that in the immediate neighbourhood of the intersection 
point, and here the term ‘continuous’ would not be appropriate. Rather, 
the width of this readjustment zone is observed to decrease, and the severity 


* In a private communication, D. R. White has emphasized another requirement 
tor the onset of Mach reflection, namely, that the stem shock must be strong enough 
to keep up with the incident shock: i.e. velocity of stem shock = (velocity of incident 
shock) * sin a,. ‘This, in turn, implies a certain pressure behind the stem, and hence 
defines the strength of the reflected shock as a function of a, if the assumption of 
pressure continuity in the vicinity of the intersection point is valid. 
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of the density gradient in it to increase rapidly near the intersection point, 
apparently condensing into a discontinuous variation right at the vertex 
of the angle between the shock and the boundary. Indeed, on the basis of 
the proposed model, a discontinuity of this sort must be. expected here, 
since, right at the boundary, the additional flow deflection must take place 
immediately. Such a variation would invalidate assumption (4) of the 
regular reflection theory; i.e. the changes in the flow in the angular region 
behind the reflection would not be negligible in comparison with those 
sustained through the shocks, regardless of how small a region about the 
intersection point is considered. 

Sharply localized subsonic variations such as this can be found in other, 
more familiar aerodynamic situations. The same type of discontinuity 
occurs in the corresponding region of the supersonic aerofoil problem 
discussed above, for example, and also in the basic problems of subsonic flow 
over a sharp convex or concave corner (see figure 12). While supersonic flows 


Figure 12. Subsonic flow over sharp corners. ‘The abrupt deflections of the ideal- 
fluid streamlines at the corner are replaced, in real fluids, by regions of 
turbulence and separation. 


can negotiate such abrupt changes in direction neatly, via oblique shocks and 
Prandtl-Meyer waves, theoretical calculations in the subsonic case become 
ambiguous right at the discontinuities in the boundaries, and infinite or 
zero Velocities appear in the results. Physically, of course, such singularities 
in subsonic flows do not occur, but are replaced by regions of separation, 
eddy formation, or general turbulence, which reflect the non-ideal nature 
of the fluid involved, and which cannot be discussed on a non-viscous 
basis. However, these zones are normally quite small, and idealized theory 
can provide an adequate description of all the remaining flow field without 
any detailed consideration of the singular region. 

Presumably, a similar dissipative process takes place in the close vicinity 
of the intersection point in the reflection problem, for here too there is 
evidence of a large velocity gradient. It would be highly interesting to 
observe directly the details of this process; but to date the experimental 
technique has not provided unambiguous description of the internal 
structure of this zone. ‘The interferograms present information on the 
integrated density along the light path, and in this sense average out any 
three-dimensional irregularities or turbulence that might be taking place 
here. Only two-dimensional phenomena such as a cylindrical region of 
separation, or an eddy, could be recognized from the interferograms; and 
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such effects, if they exist, apparently are confined to too small a neighbour- 
hood to be distinguished from the shock intersections at the boundary. 
However, several uncorrelated cases of abnormal density variations in this 
region have been observed (see figure 13, plate 4). Such effects do not 
appear to be quantitatively reproducible, indicating, perhaps, that this is 
a somewhat unstable or turbulent region. 

This inability to specify the microscopic nature of the flow-adjustment 
process does not qualify fundamentally the suggestion posed by interfero- 
grams like figures 7 (b) (plate 1) and 9 (4) (plate 2) regarding a theoretical 
resolution of the difficulty. ‘The proposal is simply that an ‘deal-fluid theory 
of reflection can be reconciled with experiment in the transition zone 
%, <% <a by the admission of a familiar type of subsonic singularity 
behind the reflected shock. 


6. SUBSONIC SINGULARITIES IN MACH REFLECTION 


The second difficulty in the reflection problem—the failure of three- 
shock theory in application to weak-shock Mach configurations—can be 
approached from essentially the same point of view as the first; i.e. a search 
is made for experimental indication that one of the basic assumptions of 
the local theory is not realistic. Again, the refraction problem closely 
parallels the situation of interest; in fact, there is really no formal dis- 
tinction between the three-shock Mach intersection, and the very special 
refraction problem of a shock incident obliquely on an interface separating 
air/air. In each case the assumed configuration and the required boundary 
conditions are the same (i.e. three shocks, intersecting at a point, out from 
which must come two parallel flows at the same pressure), and in each case, 
for x <4,, two solutions appear, one of which is discarded as physically 
unreasonable. ‘There is the distinction that in the air/air refraction case, 
one accepts the solution having a reflected wave of zero strength; for the 
Mach pattern, this branch is considered trivial. Nevertheless, the problems 
are sufficiently similar to suggest that if the general refraction pattern 
encounters flow-deflection difficulties beyond a certain angle of incidence, 
the Mach configuration may likewise suffer the same fate. ‘lhe further 
inference then is that the mechanism of adjustment should also be similar, 
namely, a process like that outlined in §5. 

Closer examination of the three-shock Mach intersection reveals certain 
differences between the flow-deflection situation here, and that in the 
reflection or refraction problems. In the latter, we needed additional flow 
deflection away from the boundary. Here we find that the flows coming 
through the incident and Mach shocks are too divergent to be compensated 
by any moderately strong reflected shock, satisfying the outlet pressure 
conditions, and we look for a mechanism to complete the deflection toward 
the slipstream. (For weak incident shocks, the solutions of the three-shock 
theory place the reflected wave at an angle (w’) greater than 90° to the line 
connecting the corner and the triple point. In extreme cases, it is placed 
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at an angle greater than 90 to the wall itself. To the present author, such 
configurations appear physically unreasonable, even in the local limit, for 
a pseudo-stationary process beginning at the corner.) 

In these cases, the interferograms (see figures 7 (c), plate 1, and 9 (c), plate 2) 
show a subsonic compression following the reflection, and simultaneously a 
subsonic rarefaction following the Mach stem, both presumably helping 
in the flow and pressure readjustment process, and both presumably 
coalescing toward discontinuities right at the triple point. Since the 
three-shock theory applied locally to the Mach intersection does not admit 
such singularities, it would then follow that the observed reflected wave 
fronts would be in disagreement with the theory. 

In committing ourselves to subsonic variations to reconcile the flow 
deflections in either the regular reflection or Mach interactions, we 
admittedly introduce certain other complications that need to be examined. 
‘These arise trom the inherent non-uniformity of any finite transient subsonic 
field, and concern (a) the pseudo-stationary nature of the configurations 
and (/) the curvature of the streamlines of the dow. In the former matter, 
Bleakney & ‘Taub (1949) have pointed out that for angles of incidence greater 
than x,, another assumption of the theory—that the process is stationary 
in the coordinate system in which the intersection point is at rest—may no 
longer be justified, and the flow may be transient. If so, this would be 
suthcient grounds for discarding the theory in this range, and the observed 
discrepancies would no longer be a problem. However, on this point the 
accumulation of experimental data is most convincing that the transition 
patterns, and even the Mach patterns that occur at yet larger angles, are 
essentially pseudo-stationary interactions. ‘The angles and shock strengths 
of the configurations appear to be independent of the absolute time, and 
the various field regions, including those that are subsonic, scale according 
vt, y ¢ within the limits of observation. In view of this, and in view of the 
added experimental fact that the dissipative mechanism, whatever it may 
be, associated with the point of singularity in a real fluid is not sufficiently 
widespread to be observed at any reasonable time after inception, it seems 
more appropriate to maintain a pseudo-stationary approach when introducing 
the suggested alterations to the theoretical model. 

The second complication involved in a subsonic region behind the 
reflected shock follows from the work of Taub (1953), who pointed out 
that if sonic signals from the corner can reach the intersection point, we 
must expect the reflection to be curved over its entire length. In general, 
the curvature of this shock will not be uniform, in which case the streamlines 
of the flow leaving it will also be curved, which clearly is at variance with the 
boundary condition for a straight wall. In this connection, the argument 
in favour of a subsonic singularity gains new strength. By analogies similar 
to those introduced above it could be shown that there is adequate precedent 
for imbuing this singularity with the property not only of deflecting the 
flow abruptly, but of straightening it as well. Again, the details of the 
‘uncurving’ mechanism could presumably be inferred from the form of 
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the density variation behind the shock, and the manner in which that 
variation condenses as it approaches the intersection point. 


7. CONCLUSION 

The present theories of shock reflection processes fail to account for 
(1), the persistence of a regular reflection process beyond a theoretical 
limiting angle of incidence z,, and (I1), the observed strength of the reflected 
shock in certain Mach reflection configurations. Direct reflection experi- 
ments are at a disadvantage in attempts to resolve these difficulties; but 
certain experiments on the shock-refraction process, the theoretical analysis 
of which is formally similar to the reflection problem, display, in corre- 
sponding regions of difficulty, pronounced subsonic variations behind the 
reflected shock. In the interval x, <« <a, it is observed that the shock 
is much stronger than might be expected near the gas interface, and is 
immediately followed by a steep, subsonic rarefaction there, whose function 
appears to be to complete the flow deflections, to apply the proper curvature 
to the streamlines, and, in some cases, to adjust the exit pressures. The 
severity and localization of this rarefaction zone increases sharply near the 
interface. In fact, the variation appears to coalesce toward a discontinuity 
at the intersection point, thereby suggesting a singularity in the theoretical 
problem. 

The close theoretical correspondence of these refraction problems to 
the transition reflection («, <« < %) and Mach-reflection situations, along 
with the observed similarities in the experimental observations, suggest 
resolution of the discrepancies in the latter by introduction of the same type 
of subsonic singularities. 

It is expected that in any real fluid such singularities degenerate into 
some sort of dissipative process, such as eddies, separation, or general 
turbulence; but no positive identification of such effects on the photographs 
has been possible. Presumably, then, these processes are confined to very 
small regions; and hence their internal details should not be significant in 
the formulation of the problem. Rather, the description of the singularities 
could be based on observations, either theoretical or experimental, of the 
subsonic regions of adjustment at finite distances from the interaction, and 
the manner in which they condense as they approach the intersection point. 


The author would like to express his appreciation to the many people 
who were kind enough to read the draft of this paper and discuss the problem 
with him: Professors Bleakney, Curtis, Emrich, Krieth, Myers and ‘Taub, 
and Dr D. R. White. He would also like to emphasize that most of the 
introductory discussion of shock reflection is by no means original, but is 
rather a digest of the experiments and papers of séveral earlier workers 
cited in the references, notably W. Bleakney, C. Fletcher, L. G. Smith, 
A. H. Taub, and D. R. White. The refraction experiments on which this 
paper is based were performed at Princeton University, and were supported 
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Robert G. Jahn, Transition processes in shock wave interactions, Plate I. 


Figure 7. Shock refraction interferograms: configurations established by plane 
shock waves incident on interfaces between air ,and methane at angles 
(a)a< a; (b) a, << a< a; (€) a >a. (Symbols are same as in figures 5 


and 6.) 
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(a) 


Figure 9. Shock reflection interferograms : configurations established by plane shock 
waves incident on a solid boundary at angles (a) a< a, 3 (b) a, << a< 4} 
(c) xX > aX). (Symbols are same as in figures 1 and 3.) 
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(a) 


Figure 11. Diamond aerofoil at three-flight conditions : (a) straight attached bow 
shocks; (4) curved attached bow shocks; (c) detached bow shock. 
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Robert G. Jahn, Transition processes in shock wave interactions, Plate 4. 


Figure 13. Isolated examples of severe density gradients behind the reflected shock. 
Three interferograms and one schlieren photograph of transition reflection 


processes. 
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The diffusion of smoke from a continuous elevated 
point-source into a turbulent atmosphere 


By F. B. SMITH 


Department of Mathematics, University of Manchester 
(Received 16 August 1956) 


SUMMARY 


The problem is to determine the downwind concentration of 
non-buoyant smoke emitted from a continuous elevated point- 
source in a turbulent airstream. ‘The velocity and eddy diffusivity 
coefficients are represented by related powers of the height above 
ground, and are independent of the position of the source. Exact 
solutions are obtained for the zero and second moments of the 
concentration distribution along lines lying in the cross-wind 
direction at ground level. In special cases, these moments may 
also be determined along lines at general height. In one such 
case the concentration is determined exactly (rather than just the 
two moments) and it is found that the cross-wind distribution 
always has a Gaussian form. If it is assumed that in all cases the 
cross-wind profile is Gaussian, a formulation for the concentration 
can be given purely in terms of the known zero and second moments. 
When the source is moving with constant velocity across the wind, 
the first moment as well as the zero and second moments is exactly 
determined, and under a similar assumption a formula for the 
concentration is found. 


1. INTRODUCTION 


A new approach is used in tackling the problem of the determination 
of the concentration of non-buoyant smoke downstream of grounded and 
elevated point-sources in a statistically-steady turbulent airstream. 

When the mean velocity is represented by a power of the height, 

u = u(z+h), (1.1) 
(z is the height relative to the source height 4; see figure 1) and the eddy 
diffusivity coefficient K is related to the velocity by the usual conjugate 
power relation based on Reynolds analogy, an exact solution for the 
concentration due to an elevated line-source is found for general values 
of the velocity profile parameter x. ‘This solution, together with a further 
exact solution for the ‘spread’ of the plume from a point-source, again 
for general x, enables an expression for the ground level concentration, 
due to an elevated point-source, to be found on the basis of the single 
assumption that this concentration has a Gaussian profile in the horizontal 
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cross-wind y-direction. ‘The restriction, that this expression is for the 
concentration at ground level only, rather than at all heights, is removed 
when the source itself is at ground level, no matter what value x may take, 
and also when the source is elevated provided that « is zero or one-half. 

The nature of these solutions and their success throws doubt on the 
theoretical basis of some previous solutions and on their range of 


applicability. 

z 
] 

u(z 

S 


Figure 1. The system of axes in relation to the position of the source, elevation h, 
and the unidirectional velocity field u(z). 


Before the mathematical analysis is given, it is desirable to elucidate 
and to add to the above statements. In particular, it will be necessary to 
discuss the nature of the eddy diffusivity coefficients, as some forms 
suggested fairly recently by other workers are fundamentally different in 
concept. How these other forms arose and developed is perhaps best 
seen in their historical setting. 

The equation of diffusion is 


coc) f,, OC 


where x is the direction of the mean stream and C stands for the concen- 
tration. ‘The diffusion in the direction of the mean stream is neglected, 
on the usual ‘boundary-layer’ approximations, by comparison with the 
other terms in the equation. It remains to specify the form of u(z), K, 
and K,,. Following Calder and others, the velocity is given by (1.1), which 
for problems of this type is a sufficiently good representation and is 
preferable on grounds of mathematical simplicity to more sophisticated 
profiles, such as the logarithmic profile based on Nikuradse’s relation, or 
Deacon’s law. The parameter « in equation (1.1) is dependent on the 
stability of the air: in neutral stability it is normally about 4, and it 
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increases as the stability increases. As we have stated, the form of the 
diffusivity coefficients is in dispute but it is true that there has been fairly 
general support for supposing that K., the coefficient representing the 
vertical diffusion process, satisfies the Reynolds analogy; that is, K, is 
proportional to the coefficient of vertical momentum transfer and thus 
(under conditions of zero horizontal pressure gradient) satisfies the 
relationship that states that the shearing stress is constant with height: 


K,= = constant. (1.3) 
Thus if wu satisfies (1.1) then K, also obeys a simple power law 
K, = K,(z+h)™. (1.4) 

At about 100 metres, the turbulence is almost perfectly isotropic; and 
even below this level the degree of isotropy seems sufficiently high to 
warrant letting K,, as Sutton (1953) has suggested, vary in an identical 
way with K,, namely, 

K, = K,(z+h). (1.5) 

This assumption may be tested by the validity of the solutions derived 
with it. Under certain circumstances it is conceivable that a more general 
form for K,, would be appropriate. By taking 

K, = K,(z+h)*"" (1.6) 
the solutions may still be found for sources at ground level and also for 
elevated sources when only the ground level concentrations are required. 
Thus, although we expect to be zero, its presence in the solutions will 
serve to indicate the influence of K,. 

The physical significance of these forms will be discussed at the end 
of the introduction in the comparison with previous work in this field. 

In tackling the three-dimensional point-source problem it is, first of 
all, desirable to try to find any exact solutions of equations (1.1), (1.2), 
(1.4) and (1.5). Section 3 deals with the exact solution for the elevated 
point-source when « = }. In this particular case it is possible to separate 
(in the mathematical sense) the variables y and z in the diffusion equation 
to give two separate differential equations, one with independent variables 
x and z, the other with y and x. The first equation is in fact the equation 
for a line source at the same x and z as the point-source, and the second 
equation gives the behaviour of the concentration profile with respect to y 
as Gaussian. It is therefore necessary to have the exact solution of the 
two-dimensional line-source problem. As will be seen below this solution 
is required not only for « = } but also for general values of «. 

In extending the work to general x, for which an exact solution could 
not be obtained, it seemed desirable to accumulate first as much exact 
knowledge about the point-source plume as possible. With this in mind, 
the approach consists in determining exact solutions for two functions of 
the point-source concentration, namely 


C, = | Cady (1.7) 
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and dy, (1.8) 


as functions of xand z. The advantages of forming the differential equations 
for these two functions are that they are now two-dimensional, the boundary 
conditions are known, and their solutions represent two of the most important 
features of the complete solution C. The first function Cy gives the total 
amount of concentration on any transverse line x = constant, z = constant; 
and, moreover, since integration of equation (1.2) with respect to y shows 
that Cy is the solution of the equivalent line-source equation, it is of primary 
importance to solve the line-source problem for this reason, apart from the 
fact that it is of course an interesting problem in its own right. It is worth 
while noting that the term representing the lateral diffusion, and thus K,, 
does not enter into the equation for Cy. The solution for the general 
elevated source is derived in §2 and is discussed below. 

The second function is the second moment of the transverse concen- 
tration profile and thus C,/C, is a measure of the ‘ spread’ in the y-direction. 

Now, since we noted that in the exact solution for « = } the concen- 
tration profiles had Gaussian transverse cross-sections, it is reasonable to 
use the exact solutions C, and C, to formulate for general « an expression 
for C which satisfies all the boundary conditions and which is based on the 
single hypothesis of a Gaussian transverse cross-section. ‘Thus 

C= X(x, 9, (1.9) 
where X(x, ) and f(x, x) can be expressed in terms of Cy and C, as follows 
(using (1.7) and (1.8)): 

f=2C,/C, and X = Cyv/(C,/2xC,). (1.10) 
The choice of the Gaussian form (1.9) is supported by the fact that, since 
the diffusion equation contains only the single derivative with respect to y, 
K,,C/oy’, with K,, independent of y, the diffusion equation effectively 
represents a process of simple diffusion in the lateral direction, a process 
with which Gaussian profiles are usually associated. 

Since it is the spread which essentially determines the rate of decay 
of C for large x, it is to be expected that the expression (1.9) will have the 
additional property that it has the correct asymptotic behaviour. The 
general elevated line-source problem has as a particular solution the case 
of zero source-elevation A, which was historically one of the first to be 
found. ‘This solution (5.1) which has received ample experimental verifi- 
cation is, in fact, included in the solution for general source elevation that 
is derived in §2 by the use of operational methods. ‘This latter solution 
involves a Bessel function but is nevertheless readily computed (especially 
if all that is required is the concentration Cy at ground level zs = —h). 


We have 


sh+h?)j2 1+2x 1 yl +2x 
2) = QO (sh+h?) exp| +h) uy 


(2x41) Kyx K,(2a+ 12x 


2up(sh +h?) 1-292 
| Ky(2%+ 1)?x |: (1.11) 


LID 


| 

| 

4 

1 

{ r 


Diffusion from an elevated source in a turbulent atmosphere 53 


which in the case of x = } has the shape shown in figure 2 for y = 0, 
x= —h. The curves for other values of x are similar in their general 
shape although the precise details will of course vary. It is interesting to 
note that as one travels in the direction of increasing x the jump in the 
concentration from zero to an appreciable value takes place very rapidly. 
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Figure 2. ‘The concentration at ground level, as a function of the distance downstream 


due to an elevated line source when « 3. 


z=0 
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Figure 3. Diagrammatic representation of the distribution of concentration along 
the vertical near the source. 


The solution has the same asymptotic behaviour for large x as the verified 
particular solution for h zero; for example, in conditions of neutral 
stability (« = 4) the concentration behaves like the (— %)th power of x. 
This ultimate similarity for all h is to be expected on realizing that the 
smoke distribution ‘ forgets’ its initial distribution as it is swept downstream. 

Near the source the concentration is asymmetrically distributed about 
the plane z = 0 (this is diagrammatically represented in figure 3). ‘This 
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asymmetry may be accounted for by the combined effect of the variation 
of K and u with height in the following way. As the vertical diffusive 
processes increase with height, the smoke is distributed along the vertical 
more rapidly above the source than below it, so that the concentration 
profile near x = 0 becomes relatively distended above the source. Initially, 
then, the concentration is greater at, say, = Athan at z = —A (as figure 10 
shows). ‘The velocity shear also accentuates the ‘asymmetry’ effect; if, 
for a moment, one neglects diffusion and imagines that a steady source 
injects smoke at a given rate into a steady stream then it is clear that the 
density of the smoke particles in that stream would be inversely proportional 
to its speed. Thus the effect of the velocity shear is to increase the 
concentration below the source, where the stream speed is lower, and to 
decrease it above the source, remembering however that the actual amount 
in any layer is governed largely by the diffusion process. 

The differential equation for C, is somewhat more complicated than 
that for Cy, owing to the addition of an extra term depending on the 
functional nature of A, and on C, itself. For zero source-elevation 4A, 
the method of solution is fairly straightforward and consists of expressing 
C, in the form 

(1+2x)?Ky x” 
By solving for C(y) and applying the appropriate boundary conditions, C, 
is finally expressed in terms of a confluent hypergeometric function and an 
allied function V with similarly rapid convergence properties : 


2 


C= x*C(n), 


(1.12) 


K,\(2\e>-2 Ke a (b—1)!(64 a—2)! 


b—1)! 


(a—1)! 
The behaviour of C, is chiefly governed by the term e~" for small x, and 
by the term x’ for large x. ‘Thus the general shape is as shown in figure 4. 
It is found that the solution of the differential equation is simplified when « 
takes the values 0, 4, } and 1. 

The problem of the elevated source does not meet with quite the same 
general success, although C, may be directly determined at all levels for 
a = 0 (as well as for « = 3, for which the exact solution has already been 
quoted) and at ground level for « = 4 and 1. (It seems reasonable to 
suppose that the solutions for these two cases could likewise be extended 
to all levels if so desired, but at present the actual method is not clear.) 
The difficulty for general x is not theoretical but rather a matter of being 
stranded with an unpleasant indefinite integral (6-6) to evaluate, and this 
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it is possible to do in terms of simple functions only in the above four cases. 
As in the case of the elevated line-source, the solutions are obtained by the 
use of operational methods. 

However, for most purposes the values of the concentration are required 
only at ground level, and this is certainly true in problems of atmospheric 
pollution. This greatly simplifies the problem since the concentration at 
ground level downwind of a point- or a line-source elevated at height A 
above ground is identical with the concentration at the same x and y, but 
at height /, due to an exactly similar source at ground level. ‘This property 
of reciprocity has been noticed by Bosanquet & Pearson (1936) in the 
case of the line-source when « = 0, and in the obvious case of u = constant, 


Cap 
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Figure +. Diagrammatic representation of the variation of the 
spread’ Cy, with distance downstream. 


B 


Figure 5. The two congruent paths Y and Y’. The coordinates of the points are : 


A= (0, 0, h), B= (0, 0, 0), C= (x, J; h), D= (x, J; 0). 


K = constant. They realized that some general theorem of the nature of 
the one given here might be true. This reciprocal theorem is rigorously 
proved in §5 and is a consequence of the fact that the solutions of the 
problem are Green’s functions. Physically the explanation of this theorem 
is not so straightforward; if, however, one assumes that the distribution 
of perturbation velocities at any point in the field is symmetrical, then the 
probability that a particle leaving the elevated source follows a certain 
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path ¥ to a point (x, y,0) at ground level is the same as the probability 
that a particle leaving the same source now at ground level will follow the 
congruent path Y’ to the point (x, y,/), since both particles pass through 
corresponding points and conditions, only in opposite order. Integration 
of all such possible paths recovers the theorem. Hence the concentration 
at ground level due to the elevated source may be determined from (1.12) 
by putting =A, that is 7 = uy h'***/(1 + 

‘The results are plotted for « = 0,3, }, } and 1 in §7, and a comparison 
between them is made. From these solutions, interpolation curves for the 
height and distance downstream of the maximum concentration are deduced 
and it is noted that the greatest variation of these factors occurs for small «. 

In the final section, §8, a simple extension of §4 is investigated. The 
problem is made more general by allowing the source to move (with the 
application to things like smoke-screen layers in mind) with velocity constant 
in the y-direction. Again C, and C, are found, relative to axes moving 
with the source. Furthermore, the function 


yCdy (1.13) 


is determined. ‘The function C,/C, = Y(x,z) is a measure of the lateral 
displacement of the centre of the plume due to the motion of the source. 
For large x, this displacement is proportional to x1*% “+2”, which implies 
that, for « > 0, the mean velocity of a group of smoke particles is 
tending to increase with time as the group is swept downstream due to 
the general deepening of the plume and the increase of velocity with 
height. 

‘lo conclude this introduction, we discuss the nature of the eddy 
diffusivity in the light of previous solutions and formulations. Naturally, 
the form chosen for A, has an important influence on the character of 
the solution; in particular, on the rate of decrease of concentration for 
large x. Briefly this is because, since this rate of decay is asymptotically 
the same at all levels, the overall decrease, which is governed by how 
quickly it spreads out laterally, is the greater the faster A, increases 
with height. For this reason the elementary solution corresponding to 
u = constant, AK, = K, = K = constant, 


O (v2 +(s—h)*)u (2+h)?)u 
C= - |+ exp » (1.14) 


is not applicable in the atmosphere. (In (1.14) z= 0 is ground level, 
h the source elevation.) We note that this solution gives a rate of decrease 
like x! whereas experiments (Sutton 1947) show that this index of x 
should be about — 1-76 (cf. equations (1.12), (1.11) which give an index 
— 1-67 for « = 4). 

Similar objections (Sutton 1953) have been raised against the solutions 
due to Davies (1950) for which u and K, satisfy (1.1) and (1.4) respectively 
and K,, obeys the power law relation 


K, = K,(z+h). (1.15) 
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(in Davies’s solution # = 0.) ‘This relation had apparently two advantages. 
First, on physical grounds, it seemed at first sight, due to some degree of 
anisotropy in the turbulence very near the ground, to be reasonable to 
suppose that AK, > K_ near the ground (that this idea leads to a wrong 
‘decay’ index is clearly tied up with the concept given above, that the 
decay is governed by the spreading out at consecutively higher levels, as 
x increases, where the turbulence becomes more and more isotropic). 
Secondly, it enabled the two variables y and s in the diffusion equation (1.2) 
to be separated in the mathematical sense in a way similar to the present 
solution for x=}. Sutton (1953, p. 283) has suggested that in fact a 
more realistic rate of decay will be obtained if A,, is taken equal to K, (as 
is stated in (1.5)). As we have seen this is borne out by the solutions 
obtained in the subsequent sections. 

Supported by the comparative failure of the above two solutions, 
different approaches to the problem have been developed. Notable amongst 
these are those due to Sutton (1934) and Davies (1954). Sutton adopted 
a statistical approach based on ‘laylor’s theorem for the standard deviation 
o(7’) of the distance travelled in a time 7 by an infinite succession of 
particles in a field of homogeneous turbulence. An expression for the 
concentration as a function of x was formulated in which the lateral 
dispersion was derived using ‘Taylor’s theorem; but except for this 
dependence on x and the consequential necessity of making the solution 
satisfy the continuity of matter equation, the expression is based on the 
solution (1.14) for constant velocity and K with height. In the case of a 
point-source at ground level, the expression is 


where of = 4A, = 4A, 


It seems, by an investigation due to Knighting (see Sutton 1953) on 
‘turbulent diffusion as a random walk process’, that the basic assumptions 
lead to a conception of the diffusive process not intrinsically assumed in 
Sutton’s derivation. ‘This must be due to the character of the assumptions 
themselves; for example, to basing the expression on the solution for 
constant velocity. ‘To make but one comment on this assumption, the 
interpolation curves of §7 indicate that in fact the concentration profiles 
are most sensitive to changes in « when « is small (near constant velocity). 
The result is that equation (1.16) represents a diffusion process which 
implies that the diffusion is being carried out mainly by those eddies that 
are smaller than, or of the same order of magnitude as, the ‘spreads’ a, 
and o.. 

In the work of Davies (1954) this concept is openly assumed. ‘The 
effect of the turbulence is represented by eddy diffusivity coefficients K,, K-; 
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path ¥Y to a point (x, y,0) at ground level is the same as the probability 
that a particle leaving the same source now at ground level will follow the 
congruent path % to the point (x, y,/), since both particles pass through 
corresponding points and conditions, only in opposite order. Integration 
of all such possible paths recovers the theorem. Hence the concentration 
at ground level due to the elevated source may be determined from (1.12) 
by putting z =A, that is 7 = u,h'***/(14 

The results are plotted for x = 0, 3, }, } and 1 in §7, and a comparison 
between them is made. From these solutions, interpolation curves for the 
height and distance downstream of the maximum concentration are deduced 
and it is noted that the greatest variation of these factors occurs for small «. 

In the final section, §8, a simple extension of §4 is investigated. ‘lhe 
problem is made more general by allowing the source to move (with the 
application to things like smoke-screen layers in mind) with velocity constant 
in the y-direction. Again C, and C, are found, relative to axes moving 
with the source. Furthermore, the function 


C,=| yCdy (1.13) 


is determined. ‘The function C,/C, = Y(x, 2%) is a measure of the lateral 
displacement of the centre of the plume due to the motion of the source. 
For large x, this displacement is proportional to x4*~%%*?”, which implies 
that, for « > 0, the mean velocity of a group of smoke particles is 
tending to increase with time as the group is swept downstream due to 
the general deepening of the plume and the increase of velocity with 
height. 

‘Yo conclude this introduction, we discuss the nature of the eddy 
diffusivity in the light of previous solutions and formulations. Naturally, 
the form chosen for A, has an important influence on the character of 
the solution; in particular, on the rate of decrease of concentration for 
large x. Briefly this is because, since this rate of decay is asymptotically 
the same at all levels, the overall decrease, which is governed by how 
quickly it spreads out laterally, is the greater the faster A,, increases 
with height. For this reason the elementary solution corresponding to 
u = constant, K, = K, = K = constant, 


(y?+(s—h))u (v?+(s+h)*)u 


is not applicable in the atmosphere. (In (1.14) s =0 is ground level, 
h the source elevation.) We note that this solution gives a rate of decrease 
like «~! whereas experiments (Sutton 1947) show that this index of x 
should be about — 1-76 (cf. equations (1.12), (1.11) which give an index 
— 1-67 for « = 4). 

Similar objections (Sutton 1953) have been raised against the solutions 
due to Davies (1950) for which u and K. satisfy (1.1) and (1.4) respectively 
and K,, obeys the power law relation 


K, = K,(z+h). (1.15) 
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(In Davies’s solution # = 0.) This relation had apparently two advantages. 
First, on physical grounds, it seemed at first sight, due to some degree of 
anisotropy in the turbulence very near the ground, to be reasonable to 
suppose that A, > AK. near the ground (that this idea leads to a wrong 
‘decay’ index is clearly tied up with the concept given above, that the 
decay is governed by the spreading out at consecutively higher levels, as 
x increases, where the turbulence becomes more and more isotropic). 
Secondly, it enabled the two variables y and z in the diffusion equation (1.2) 
to be separated in the mathematical sense in a way similar to the present 
solution for x= 3. Sutton (1953, p. 283) has suggested that in fact a 
more realistic rate of decay will be obtained if A, is taken equal to K, (as 
is stated in (1.5)). As we have seen this is borne out by the solutions 
obtained in the subsequent sections. 

Supported by the comparative failure of the above two solutions, 
ditferent approaches to the problem have been developed. Notable amongst 
these are those due to Sutton (1934) and Davies (1954). Sutton adopted 
a statistical approach based on ‘Taylor’s theorem for the standard deviation 
o(7) of the distance travelled in a time 7 by an infinite succession of 
particles in a field of homogeneous turbulence. An expression for the 
concentration as a function of x was formulated in which the lateral 
dispersion was derived using ‘Taylor’s theorem; but except for this 
dependence on x and the consequential necessity of making the solution 
satisfy the continuity of matter equation, the expression is based on the 
solution (1.14) for constant velocity and A with height. In the case of a 
point-source at ground level, the expression is 


20 
C(x, y, = 7A, A. t ’ (1.16) 


where of = $A, of = 44, 


It seems, by an investigation due to Knighting (see Sutton 1953) on 
‘turbulent diffusion as a random walk process’, that the basic assumptions 
lead to a conception of the diffusive process not intrinsically assumed in 
Sutton’s derivation. ‘his must be due to the character of the assumptions 
themselves; for example, to basing the expression on the solution for 
constant velocity. ‘To make but one comment on this assumption, the 
interpolation curves of $7 indicate that in fact the concentration profiles 
are most sensitive to changes in x when « is small (near constant velocity). 
The result is that equation (1.16) represents a diffusion process which 
implies that the diffusion is being carried out mainly by those eddies that 
are smaller than, or of the same order of magnitude as, the ‘spreads’ a, 
and o.. 

In the work of Davies (1954) this concept is openly assumed. ‘The 
effect of the turbulence is represented by eddy diffusivity coefficients K,,, K-.; 
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for a source at ground level, the form of A. remains as in equation (1.3), 
whereas K, is made a function not only of z but also of y, namely 


K,, = K,2%y'*, (1.17) 


so as to make the diffusion more rapid as the plume spreads laterally. 

‘This form of A,, retains both the advantages of K,, in equation (1.15), 
and moreover the solutions of (1.2) have a more closely correct rate of 
decay for large v; the form is nevertheless open to considerable theoretical 
doubts. (Assumably, for an elevated source, AK. would likewise depend 
on the coordinates of the source.) 

In considering this fundamental assumption, a clear distinction has to 
be made between two different kinds of problem, namely ‘one-particle’ 
problems and ‘two-particle’ problems. ‘Thus our problem, in which we 
require the expectation concentration downwind of a source, that is the 
probability that a single particle leaving the source will arrive at the point 
in question, falls into the first category. On the other hand if we are interested 
in, say, the expectation of the maximum density in the plume at any time 
some fixed distance d downstream, regardless of where this maximum may 
be in the plane x = d, then we clearly have to know the way in which one 
particle is tending to diffuse away from its neighbour; this problem falls 
into the second category. Now one-particle problems of this kind are 
problems in which the particle is subjected to displacements by all eddies, 
large or small, whereas in two-particle problems those eddies whose 
characteristic size is greater than the separation of the particles are not of 
primary interest in the problem because they merely shift the two particles 
equally and do not contribute to the factor required, namely their separation. 
And so, if we consider the point-source plume at any particular time, those 
eddies that are smaller than the width of the plume are spreading the plume 
out relative to its own axis whereas the larger eddies are bodily shifting the 
plume from side to side and up and down, and are thereby introducing 
‘intermittancy’ in the catch at any recording instrument, similar to that 
in the turbulence at the edge of a wake. 

Thus in two-particle problems it is indeed necessary to choose eddy 
diffusivity coefficients dependent on the position of the source. But since 
the present problem is a one-particle problem in which all the eddies are 
everywhere contributing to the diffusion process, that is to the mean catch, 
it is necessary to choose coefficients independent of the source position, as 
in fact has been done in this paper. 


2. ELEVATED LINE-SOURCE SOLUTION C, 


The infinitely long line-source lies across the wind in the horizontal 
y-direction (figure 1) and is elevated a height A above ground level. When 
z is the vertical coordinate measured with the source as origin the velocity 
and the eddy dittusivity coefficient are taken to satisfy equations (1.1) and 


\ 
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(1.4) respectively. ‘The expectation concentration then satisfies the equation 
of ditfusion 
+h) = | | (2.1) 
with the following boundary conditions. 
(i) In the plane of the source x = 0), C is zero except at the source, | 
where it is infinite: in ‘delta function’ notation, C = O 6(2z). 
(ii) ‘The ground is impervious to smoke: on z = —h,K.eC/dz = 0. 
(iii) ‘he concentration dies away at great heights: C > Oasz-> &. 
(iv) ‘The flux across any plane x = constant is independent of the } (2.2) 
value of x, since no smoke is lost or created for x > 0: 


| Cudz=Q. 


—h 
(‘The last condition is not independent of the others.) 


For the sake of simplicity the mathematical manipulation is carried out 
with coordinates related to the real ones as follows: 


hay, (2.3) 


Uy 
Ky 
(the real coordinates are placed first). ‘lhis results in elimination of the 
constants uy, Ky and h from (2.1) and (2.2). 
The method of attack is to use Heaviside p-operators whereby the 
ditferential operator ¢/ex is replaced by p; and C(x,z) is replaced by 
C(p, 2), where 


z— zh, 


C(x, z) = C(p,3) dp. (2.4) 
Equation (2.1) becomes 
las = —(2+ 1)*"p0 A(z). (2.5) 
ovr 


The right-hand side may be put equal to zero provided the equation is 
solved separately in the two regions zs > 0 and z < 0. The equation is 
then Bessel’s equation, the solutions of which, satisfying (2.21, ii, ili), are 

C(p, 3) = > 9), \ 


C(p,2) = (-1<2z<0)f 9 


where r= + 1) +2002 (2.7) 


(throughout this paper K,(z) is defined as }zcosecv7{l_,(z)—J,(2)}). 
Also, the two solutions must be equal on z = 0; the gradients do not, 


however, have to be continuous. ‘Thus if r = q on z = 0 we have 
A(P) = (2.8) 
Finally, applying (2.2iv) and using the recurrence relation (Watson 1944, 
p. 80, eqn. (20)), 
2vp 1/(1+2a) 
D(p) = ) Ka (2.9) 


1+ 2x 


| 
3 
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and similarly for A(p). ‘The operational forms (2.6) have to be interpreted 
to recover C(x,z) by application of the integral (2.4). Clearly the two 
forms are very similar and, in fact, they have identical interpretations. 
Thus, considering only the solution for z > 0, 


2 2 


2 
x K avp dp (2.10) 
where a = (z+ 1)''°?”*. ‘This may be evaluated by deforming the contour 
to one whose significant part is a loop round the negative real axis (figure 6). 
Since the integrand has zero residue at p = 0 and gives no contribution 


C +100 


p- plane 


C-j00 


Figure 6. ‘The contour from C—-i% to C 1% is deformed into the semicircle at 
infinity and a loop round the negative real axis and the origin. 


on the circle at infinity, (2.10) may be represented as the sum of two 
integrals along the negative real axis, which may be combined to give 
: 1+2« 
(x, 3) = 9 O(z+ | I 2a(4q) x 
0 
x (2.11) 
This integral is a special case of Weber’s second exponential integral 
(Watson 1944, p. 395) 
. 1 a* + b ab 
*xp(— p???). = exp{ 2. 
exp( — p?t?).J (at).J ,(bt)t dt ( (=): (2.12) 
which, when applied to (2.11), and reverting to the real variables (2.3), 
yields the final solution 


(1422) Kgx 


K,(22 + 1)?x 


x K,(2%+ 1)2x | (2.13) 


C 
(x, 2) 


Diffusion from an elevated source in a turbulent atmosphere 61 


The general properties of this solution have been discussed in the 
introduction. 
3. ELEVATED POINT SOURCE « = } 

This is the single exact solution for the complete point-source problem. 
The particular virtue about « = }, » = 0 is that u and K have the same 
functional form. ‘This enables the variables y and z to be separated in 
the diffusion equation. 

The equation of diffusion is (1.2) with u, A, and K, given by (1.1), 
(1.4) and (1.5), respectively, with boundary conditions identical with (2.2) 
except for (iv) which is changed to 


| Cudyds = 0. (3.1) 


As in §2 it is found convenient to use related coordinates : 


x, >hz, y O > Q (3.2) 


(the real coordinates are placed first), where K,, as used in (1.6), is equal 
to K, when » = 0. The equation of diffusion is thus 


ac aC 
(s+ 1)! | +(2+ 1)! (3.3) 
A solution of the form 
C = X(x,y) Y(n), (3.4) 


is sought, whereby (3.3) can be separated into two distinct equations linked 
only by the parameter m: 
1 OX m 


ey ( =) 
3.6 


The solution satisfying the boundary conditions is found by putting m = 0. 
Then (3.6) can be solved to give 

oY 

on 
and since by inspection of (3.3) we see that 0C/dy is also a solution then 


so iS 


1 
X(x, 2) eu (3.7) 


where X(x,2) satisfies (3.5) with m=0. ‘This is exactly the same 
differential equation as (2.1) with a=}; also X(x,z) must satisfy the 
same boundary conditions (2.2) except that the flux QO must be altered 
to O/2v7. ‘Thus the solution, reverting to the real variables (3.2), is 

O (sh+h*)!4 y?+h?+(2+h)? sh+h? 

(Kyx)3? 


x/Uy 


ss 
ae 
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In the generalized case when K,, satisfies (1.6), the exact solution is found 

in an exactly similar manner for the case 
a= ,). (3.9) 

The solution is then 

O (sh+h*)j? 
2(1+2%)4/(7K,) 
1(1 + 2x)?(Ko/K,)y? + (2+ + hi +e 
_ 4 0 
exp| (1+ 22)?Ky x/ty 


C = 


X (1+ 22)?K 3. ) 
O° 
The important thing to note in this section is the Gaussian form of Y, as 
. . P . 
this is the basis of the subsequent sections. 


4. POINT-SOURCE AT GROUND LEVEL 
With the variables changed, for the sake of simplicity, as follows, 


the equation of diffusion (1.2) is 
ec 0 ac 
The approach of this section has been discussed in the Introduction and 
involves determining exact solutions for Cy (see (1.7)) and C, (see (1.8)) 
which may then be utilized to find the concentration C by the formulae 
(1.9) and (1.10). 


C,, being the solution of the line-source problem, is already determined 
($2), and in the special case of zero source-elevation h is 


where O 
(1422)! — /(1+2x)]! (+4) 


The second function, C,, satisfies the equation 


ac; g! 20 
Ox cs (1+ 


Writing 


7 1+ 2x)?x 1+ 2a 1+2x’ 
then the solution of (4.5) is obtained by expressing C, in the form 
Cy = 2A(1 + 2a)? (4.7) 


where G() satisfies the ordinary differential equation 


2G 


x 
2 
| 
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The homogeneous equation is the equation for the confluent hyper- 
geometric function and has the two independent solutions (excluding 
a = 0) 


—a)! (a—1)! 


We choose U(b;a;n) for the second solution since it is o(e”) at 7 = 2 
(Jeffreys & Jeffreys 1946, p. 579). ‘The particular integral of (4.8) can be 
expressed, by the usual formula given in the theory of ‘variation of 
parameters’, in terms of the right-hand side and the Wronskian 


a-—1)! 
so that the full solution is 
(b—1)! 
C, = + x 


x {U(b;a;7)[1, + (4.11) 


where J, and /, are indefinite integrals arising from the particular integral, 


= | do, 
I, = | [ ») dn, | 
“0 
and A and B are constants to be determined. 
The boundary conditions that (4.11) must satisfy are (2.2i, ii, iii). 
Thus C, 0 as 7 which in terms of G(y) implies that G = 0(e”) 
for large y. In this limit, the second term is dominant and would be O(e”) 


(4.12) 


unless B= —I,(«). I[,(«) may be evaluated by use of the integral 
representation of U(b;a;n), namely 
U(b;a;n) = |, + dv, (4.13) 
to give (6—1)!(a+hb—2)! 
I(o)= -—B= (4.14) 
Applying the zero-flux condition on z = 0, 7 = 0, which is 
ec 
on 0, (4.15) 


to (4.11), the only non-zero contribution is the term containing the factor 


= U(b;a; 7). 


And so, to satisfy (4.15), we have 
A=0. (4.16) 


| 

4 

| | 

| 
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Hence 
b—1)! ; 
| Usain) te" 9) dy + 
+ |, ni t-2e-" U(b; an}. (4.17) 


This is further simplified by expanding the integrands in powers of », 
integrating term by term, and then collecting like powers of 7 to give the 
rapidly convergent series 


(b—1)!(6+a—2)![(b-1)! 
where = (2b+r—1)! 


V(b;a;n) = > 


The evaluation of this solution is sufficiently simple to warrant its use. 
Generally no more than eight terms are required to give the solution correct 
to four places. 
In terms of the real variables (4.1), the ‘spread’ function is 


= OC (1 4. yy 
C, 20 MR) 20.) (a 1)!(25 i 1)! Xx 


(6-1)! 


? 
Ug 2 l+a 
(1+22)?K,x 1 + 2 +2 
Four special cases are known in which the solutions of (4.5) can be obtained 
very much more directly. ‘These solutions are obtained in terms of simple 
functions and can be shown to be included in the general form (4.20). 
‘The four cases correspond to values of x = 0, }, | and 1, with p = 0: 


OK, 


Uys 


J 
K 
Cy = 2 (4.24) 
0 
(this being the exact solution of §3); 
20 1 
n 
I*(y,€)=1-— | | dy (4.26) 


0 


(Pearson 1922), 


bi 


ec 


2 /3\15 x25 
(4.23) 
4s (3)! Duy = 
th 
| tw 


) 
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We noted in the Introduction that the rate of decrease of concentration 
as x increases for large x is considered to be a valid test of any proposed 
expression. Equations (1.10), (4.3) and (4.20) show that 


Com — (1+2a) (4.27) 
so when » = 0 and « = } (the expected value in neutral stability) 
(4.28) 


which agrees with those values which have been obtained by experiment 
(Sutton 1947). 


5. RECIPROCAL THEOREM 


The reciprocal theorem, which is of vital importance in extending the 
work to the problem of the elevated source, may be stated as follows. 


The concentration at x’ due to a source at x”, with the flow in the 

positive x,-direction, is equal to the concentration at x” due to an 

identical source at x’ when the direction of the flow is reversed. 
The proof of the theorem runs thus. 

Consider the partia! differential equation 


d/,- 9G 0G 


(with summation over repeated suffixes), where G is a Green’s function 
with centre x” and which satisfies - boundary conditions 


OX; 


= 0 
aG on the ground Wy 
OX; J 


(n; are the direction cosines of the normal to the bounding surface—for the 
horizontal ground ; = (0,0,1)). Let this equation be represented by 
LG = 0 where L is an operator. Also let the adjoint equation, with the 
same boundary conditions, 

K,,=—) — —(B,G,) =| 5. 

(B; G,) = 0, (5.3) 
be represented by MG,=0. G 


function of the adjoint 
equation with centre x’. Then 


| dx 
Jy 


Ox. 


0 


| I. GG,, B,n, aS, (5.4) 


by the of Green's and the divergence theorem. Let 
the bounding surface S be the ground, the hemisphere at infinity and the 
two small spheres excluding the singularities x’ and x”, Since B;n; = 0 
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on the ground and G~1~,r (r= distance from source x”) the second 
integral makes no contribution. 

The first integral makes a contribution only on the small spheres: 
the first term on the one around x”, the second on the one around x’. 
Provided these spheres are sufficiently small, G,, is constant on the one 
at x” and G on the one at x’. Since the flux from the two sources are 
equal, the right-hand side of (5.4) reduces to 


O[G,(x", x’) — G(x’, x”)], (5.5) 
where O is this flux. ‘Thus 
G(x’, x") = G,(x’,x’); (5.6) 


in words, the Green’s functions of the equation and its adjoint are the 
same if the variables are interchanged. 


Now if 
0) 0 
K,(z) 0 | and B; = (u(s),9,0), (5.7) 
0 O | 


then LG = 0 is the diffusion equation with flow in the positive direction 
of x and MG, =0 is the diffusion equation with flow in the negative 
direction of x. Also, G(x’,x”) represents the concentration at x’ due to 
a source at x” with flow in the positive direction, and G(x", x’) represents 
the concentration at x” due to a source at x’ with flow in the reversed 
direction. We have G(x’,x”) = G,(x’,x’) which is the theorem to be 
proved. 

This theorem, as discussed in the Introduction, enables the ground 
level concentrations to be determined downwind of an elevated point-source 
from the concentrations at arbitrary level downwind of a point-source on 
the ground. In particular it is used to determine C, at ground level 
(everything above applies to Cy and C, as well as to C). For source-elevation 
h this is obtained from (4.20) and (4.21) by putting s = h. 

6. ELEVATED POINT-SOURCE * SPREAD’ FUNCTION C, 

It would be desirable to find C, for all heights (although of course it 
is the ground level value in which we are particularly interested). ‘The 
complete solution has already been found in $3 for x= 5. ‘This section, 
following an approach similar to that used in §2 to find Cy, yields the 
complete solution for x =. It seems possible that the cases x = } and 
z = 1 could also be treated in this way. 

With the transformation of variables (3.2) together with changing C, 
to h®(,, the modified C, satisfies the equation 

(s+ = +3)! +2(14+s)'°C,, (6.1) 
when is zero. In p-operational form (see § 2), this becomes 
ec, 


— 
1+2 02 


d 

‘ 

4 


2) 
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since C, = 0on x = 0. As before, put 


? 
r= 14 20 + 1) +201) 2 (6.3) 
and let 
1 + 2 
20 6 
Then (6.2) becomes 
by the application of (2. 2.8) and (2.9). The complementary functions 


are the same as in §2. gh trouble arises in attempting to evaluate the 
particular integral. For z —- 0, say, the particular integral is 


0 
Al | 9 ons) ds. (6.0) 
- 


These integrals may be evaluated in terms of simple functions only for the 
four values of x: «= 0, 4, and 1 (cf. §4). The solution for } has 
been carried through and is in complete agreement with §3. ‘The two 
cases of x=} and x= 1 have been explored; however in taking the 
interpretations of the solutions obtained, it was found necessary to put 
> = —1 (corresponding to ground level). ‘Vhe results were in full agree- 
ment with the solutions obtained by application of the reciprocal theorem 
to (4.23) and (4.26). 

For x = 0, the analysis can be carried through without any restriction. 
Since x = 0 is of particular interest, being close t» z = 4, and also because 
it is typical (although quite ditferent in detail) of the other three solutions, 
the solution will be given below. 

When x = 0, equations (6.4) reduce to 


4(p) = 


(6.7) 
| 
A'(p) Kl2sp) 
and (6.6) reduces to 
= AK,(r) | Ky ds—Al(r) | s'K2ds. (68) 


A similar expression is obtained for s <0. If C, and @, represent 
modified Bessel functions of purely imaginary argument (/,, and A, cos nz) 
then 


re 


| FC (s) ds = [3C(r)@,(r) —2C (6.9) 
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(Watson 1944, p. 136). This relation is used to express the particular 
integrals (6.8) in terms of products of three Bessel functions. By re- 
arranging the terms, the Wronskian relationship (Watson 1944, p. 80, 
eqn. (20) 

K + L()K, alr) = Ur (6.10) 


may be used to give 


(6.11) 
I,(r) 
Thus the solution which has the correct behaviour at z = & is 
(2 >0), | 


4 I, 


Where L(p), F(p) and G(p) have to be determined from the ooniliien: 
(i) no flux across s = —1, r= 0, 1.e. 0C,/er = 0, giving F(p) = 0; 
(ii) on x = 0 the two solutions have to be equal, so that, if r = g = 3vp 
on z = 0, we have 


Cy = — }K,(r) + 
2 | 


This equation (6.13) has to be reduced to a form in which it may be 
interpreted; we are limited to reducing the two relations to sums of pairs 
of Bessel functions of equal order, or the derivatives of such pairs. ‘The 
clue, as to what form to look for, is given by the term ( )r3Q/,(q)K,(r). 
It suggests that the form should contain derivatives up to and including 
the third. By use of the ditferential equations satisfied by the Bessel 
functions and also the recurrence relations such as 


2 
I(q) = I.(q)- (q) = = (6.14) 


the equation (6.13) may be simplified to 


} 
(6.18) 


[3 


J 
The interpretation will be carried out for 3 — 0, The method is exactly 


j 
J 
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similar for z > and it can be seen that the interpretations would finally 
be identical, as they were in § 2. 
The second term has the interpretation 


Tt) op 
and on integrating by parts this is 
» *U+ta 
J, = 20(3+1)! dp. (6.17) 


Following the same process as in § 2, we find 


I, = 20(z +1)! *exp| (6.18) 


x 


The first term in (6.15) (s < 0) has the interpretation 


1 -C+iw a? ever 
which again may be integrated by parts to give 
I, = 50 — + dp (6.20) 


C+ia 


Again following the same process as in §2, we > finally obtain 


2= 39 2 


for all z. This solution reduces to (4.22) with y = u)h/K,x when the real 
variables are substituted back into (6.22): 


. OK, 
_ 3+2h (sh +h?) 
where é= p= 
and when s is put equal to its value at the ground s = —h. Equation (6.23) 


is used together with (2.13) in (1.9) and (1.10) to plot the concentration 
at three different levels (figure 10) in the next section. 


7. RESULTS 


The results of the previous sections are perhaps best presented in 
graphical form, when the differences between the curves are most readily 
appreciated. Sinke the graphs are more or less self-explanatory it will 
suffice to tabulate the main points in a very concise way: 

The profile for the elevated line-source, when x = +, has been given 
in figure 2 in §1. 
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For the point-source problem, the reciprocal theorem has been applied 
to equation (4.20) to give the ground level concentration when x = }. 
‘This profile and the four simpler cases ((4.22) to (4.26)) are given on the 
same graph (figure 7), and are sufficient to deduce, for general x, inter- 
polation curves for the main features, namely, the magnitude of the maximum 


on T T T T T 
C 
he x 
Ko 
Up Al + 2x 
Figure 7. ‘The concentration at ground level along the line 4 0, as a function of x, 


the distance downstream, due to an elevated point source. ‘The five curves 
correspond to values of a — 0, 4, 4, 4 and 1. 
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Figure 8. Interpolation curves showing the maximum concentration C,,, at ground 
level and the downstream distance x,, from the elevated source to the position 


of that maximum, plotted as functions of the parameter «. 
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concentration and its distance downstream (figure 8). ‘hese curves indicate 
that the profiles are most sensitive to changes in x when « is small. The 
combined effect due to the variation in the shear and in K as « increases 
can be summarized as follows: 

(i) the maximum goes up; 

(ii) the maximum occurs earlier ; 

(iii) the maximum is ‘ peakier’ ; 

(iv) the plume strikes the ground earlier ; 

(v) the concentration falls off for large x more slowly. 
Figure 9 shows the variation of concentration for x = 0 at three different 
levels: at ground level, at source height, and at twice the height of the 
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Figure 9. ‘The concentration on y — 0, downwind of an elevated point source 
when a 0. It is plotted at three heights: at the ground —-—, at source 
height —-— and at twice the height of the source:+++++. 


0.4 
su,(2+h) 
0 Wy h= | 


0-2 0.4 0.6 08 1.0 1.2 
Kix 
Ug ht 
Figure 10. Both curves give the concentration when the eddy diffusivity coefficient 
is constant with height. he ‘ peakier’ of the two curves applies when the 
velocity varies linearly with height, the other when the velocity is constant. 
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source, s =h. Initially, as we saw in the Introduction (figure 3), the 
concentration is greater at z =f than at z = —h. But later, due to the 
reflective nature of the ground and to the greater lateral and vertical 
diffusion above the source, the concentration at z = A starts to fall whilst 
its value at the ground continues to increase, rising to a maximum con- 
siderably in excess of that previously experienced at z=. Eventually, 
as in the case of the line-source (see § 1), the three profiles asymptotically 
approach one another for large x. 

Finally figure 10 shows the effect of the shear when the eddy diffusivity 
is kept constant. For one curve the velocity shear is constant, conforming 
with the solution « = 1, while the other curve represents the solution for 
u = constant. The effect of the shear is summarized thus: 

(i) the maximum is increased ; 

(ii) the maximum occurs earlier ; 

(iii) the maximum ts ‘ peakier’ ; 

(iv) the rate of decrease of C with x for large x is unaltered. 


8. MOVING POINT-SOURCE 
The idea behind the previous sections can be carried over and extended 
when the source is moving with constant velocity v in the transverse 
y-direction (figure 11). ‘he source may be elevated at height / or at ground- 
level, since the reciprocal theorem of §5 is still applicable. The solution 


v 
Figure 11. he source is moving with velocity v along the ground in a direction at 
right-angles to the unidirectional velocity field u(z). 


will therefore be sought for the ground-level source. It is found convenient 
to work with axes moving with the source so that the problem is that of a 
stationary point source in a wind field (u(z), v, 0). 

With u, K. and K,, following (1.1), (1.4) and (1.5) and by the trans- 
formation of variables (4.1) with K, = Ko, the diffusion equation is 


aC 
=| ] +, (8.1) 
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This equation can be used to find the exact solutions for three functions 
of the concentration. As before, we may determine Cy and C, ((1.7) and 
(1.8)), and since the distribution is no longer symmetric in the y-direction 
the third function is chosen to be 


C,= | yC dy. (8.2) 
Now if y = Y(x,z) denotes the lateral displacement of the centre of the 


plume then 


| fy— ¥(x, z)}C dy = 0. (8.3) 

yCdy=Y¥(x,2)| = Y¥(x,z)Cy. (8.4) 
Therefore Y(x, 2) = C,/C,. (8.5) 


Equation (8.4) gives the physical interpretation of the function C,. The 
function C, is still connected with the spread of the plume. 
The spread is equal to C#/Cy, where 


cr= dy. (8.6) 
Thus 
Y*(x, 2). (8.7) 


The concentration may be represented by an expression equivalent to 
(1.9): 


C = X(x, z)exp—[{y— (8.8) 
X(x, 2) = Cy J (x3). f(x,2) = (8.9) 


The first function C, is unaffected by v and therefore is given by (4.3), 
(4.4) and (4.6): 


The function C, satisfies the differential equation 


eC oC. v 

Putting 

C, = A(1 + 2a)"%e-"@(n), (8.12) 
0 

equation (8.11) becomes 

ee a 0@ a 

+ 1) —-@= 8.13 

ox? (; 7 


the solution of which is 


€(n) = Ae'+ Be" | + — nt edn. (8.14) 
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Now C, must satisfy the no-flux condition at the ground: 9“(6C,/ey) = 0. 
Thus B=0. (8.15) 
Also C, + 0 as thus 
W hereby = le dn, (8.17) 
where v A 
M= (8.18) 


( 
So that for large x, or for s = 0 (reverting to the real coordinates), 
(14+22)*? 
L 


The shape of Y(x,z) on = 0 for = 0, 4, 


Y(x, 2) (a—1)! x“. (3.19) 


!, 1 are plotted in figure 12. 
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Figure 12. he displacement V(x, 0) of the centre of the plume at ground level, 
relative to the position of the moving source. 


‘The third function is rather more involved due to the effect of the velocity z. 
However, we may write C, = C,,+Cy. where C,, is the same as (4.20), 
in the real coordinates, and does not involve 7. C,, is the perturbation 
term and is found by a method very similar to that used in §4. It satisfies 
the differential equation 


Ox ( Ca ) K, 1 ( ) 


in the variables (4.1). 
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Diffusion from an elevated source in a turbulent atmosphere 


Putting 
= 2— (8.21) 


then (8.20) becomes 


a cx 7 
t (< 2- x = —7*e" | dy. (8.22) 


ey” 1) 07) 


The solution, satisfying the zero-flux condition on 7 = 0, is 
2a—1)![,_. 

x= | | F (2a; | t*le~! - 
“0 


a; B+ | 2U(2a;a;y) te+tdtdy tI. (8.23) 
J0 


must be determined from the condition that C, +0 as By a 
similar process as in § 4, 

(2a—2)!(3a—2)! 
a(ta— 2)! 
‘The particular integral in (8.23) may be simplitied by integrating by parts. 
It is found desirable to express U(2a;a;1) in terms of the corresponding 

,F\-functions as before: 


B= 


(8.24) 


l 
1) “ F\(1+4;2-a;y) | (2a;a;) dy- 
= 0 


“0 


| le" | F (2a; a; t) dtdy 


~ 0 
—,F(2a;a;n) | | F\(14+a;2-a;t) dtdy. (8.25) 
0 
Fortunately this can be simplified a great deal. Consider the first two 
terms. Integrating term by term and summing, we find that these terms 
reduce to 
a(2a—1) 
(see (4.26) for /*). ‘The other two terms aiso simplify and can be expressed 
in terms of the function ](b;a;7) defined in (4.19). Integrating term 
by term and collecting like powers of 4, we finally obtain 


2a—1)!(3a—2)![(2a—2)! 
a(4a— 2)! | 1)! 


I*(y,a—1) (8.26) 


a(2a—4) 
from which C,. may be obtained using (8.21) and thus, with the use of 
(4.20), C,. It should be noted that Cy, is proportional to ‘This 
perturbation term increases the spread over the stationary-source case, 


I*(y,a—1), (8.27) 


and for large x this increase is proportional to x?" (also on = = 0). 
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9, CONCLUSIONS 


‘The main conclusion is that this approach has led to a fair measure 
of success in tackling the general problem of an elevated point-source. 
The most important results are represented by the equations (2.13), (3.8), 
(4.20), (5.6) and (6.23) and the graphs. Between them, they give an almost 
complete picture of the plume under varying conditions. The reciprocal 
theorem of §5 has proved highly useful, but it should be pointed out in 
warning that this theorem only holds when the form of the eddy-diffusivity 
coefficients are such as have been taken in this paper (namely, independent 
of the source position) and cannot always be applied to solutions found 
under different assumptions, such as those of Davies. However, this is 
in no way serious, since, as has been indicated in the Introduction, the 
form of A used in this paper is considered to be theoretically preferable 
in the region in which the A-theory can be applied. 
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SUMMARY 


Certain two-dimensional, laminar boundary layers are con- 
sidered whose streamlines are closed. ‘The speed at the solid 
boundary is supposed uniform, the boundary outline being 
stationary, and the speed in the boundary layer is supposed to 
differ only slightly from that of the boundary. A formal solution 
is then obtained for the motion in the boundary layer. The analysis 
confirms that a closed boundary layer may exist and yields a 
condition needed to determine the inviscid motion. ‘The condition 
is extracted in a simplified but approximate form and two examples 
of its use are given. A further class of closed boundary layers, 
namely those for which the pressure is uniform, is also considered. 
For this class the condition needed to determine the inviscid motion 
may be derived in a form both simple and exact. 


1. INTRODUCTION 


In the steady motion of a slightly viscous incompressible fluid the 
direct effect of local viscous forces may generally be neglected except in 
the neighbourhood of certain singular surfaces. ‘The indirect effect of the 
viscous forces that act in these layers is, however, often appreciable through- 
out the whole region of motion. In such cases, the motion of the fluid in the 
‘inviscid’ region cannot be determined independently of the fluid in the 
viscous shear layers. 

In two-dimensional motions, indeterminacies which arise from the 
complete neglect of viscous action are usually of one of two kinds. Either 
a shear layer is shed from a bounding surface and the inviscid motion is 
indeterminate because the position of this shear layer is not known. Or 
a shear layer, whose position is known or can be regarded as known, is 
closed and the inviscid motion is indeterminate to the extent of certain 
constants. It is with this second kind of indeterminacy that we shall be 
concerned. 

For irrotational motions, such indeterminacies generally take the form 
that the circulation round each closed boundary is unknown. For instance, 
if fluid of infinite extent streams past a circular cylinder supposedly 
rotating so rapidly that the boundary layer does not separate, then the 
inviscid motion is undertermined to the extent of its circulation about the 
cylinder. 
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For rotational motions, the vorticity may also be undetermined. As 
a typical example may be cited the motion between two rotating circular 
cylinders whose axes are set apart. Here the direct action of viscosity may 
be supposed confined to boundary layers at the cylinders. Since the 
streamlines of the inviscid motion are closed, the vorticity of the inviscid 
motion is uniform (Batchelor 1956a). ‘The value of this vorticity and the 
value of the (inviscid) circulation round some given circuit enclosing the 
inner cylinder are then undetermined. Flow patterns similar to this have 
recently been suggested as representing the motion in closed regions divided 
from a main flow by a separated boundary layer (Squire 1956, Batchelor 
1956b). 

To resolve the indeterminacies in steady motions of this kind, it is 
necessary to examine shear layers whose streamlines are closed. Only one 
solution (Squire 1956) for the velocity distribution in a closed shear layer 
has so far been reported. In this case the motion in the shear layer was 
calculated using a Rayleigh-type approximation, the convective velocity 
in the momentum equation being replaced by the tangential velocity at 
its outer edge. 

In the following a class of two-dimensional boundary layers is considered 
for which the motion in the boundary layer can, in principle, be solved 
exactly. ‘The closed solid boundary is supposed to move tangentially 
with uniform speed, and the speed in the boundary layer is assumed to differ 
only slightly from that of the boundary. A formal solution may then be 
derived by expanding the velocity in a power series in a small parameter 
representative of the small ditferences of speed through the boundary layer. 
As may be expected, the analysis of the closed boundary layer enables the 
associated indeterminacy of the inviscid motion to be resolved, in principle, 
exactly. 

‘Those closed boundary layers in which the pressure is uniform, the 
boundary speed being non-uniform, are also considered. For these boundary 
layers the information needed to determine the inviscid motion can be 
extracted without solving in detail for the motion in the boundary laver. 


2. SERIES SOLUTION 

We start by deriving a series solution for closed boundary layers of the 
first class in which the speed is almost uniform. 

‘The motion in the boundary laver is assumed steady and two-dimensional, 
the fluid being incompressible and of kinematic viscosity v. Each streamline 
is assumed to circumscribe the closed boundary and, to the boundary layer 
approximation, to have the common length 271. ‘The velocity at the 
boundary is supposed tangential and of uniform magnitude uw). For 
detiniteness, the closed boundary may be visualized as an inextensible band 


which lines a fixed cylinder and moves transversely to its generators. Let 
vL be arc-distance along the boundary from some fixed reference point. 
Further let y (vu) Lj and uy u(x,is) be the stream function and tangential 
component of velocity in the boundary layer. With x, & as coordinates, the 
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equations that govern the motion in the boundary layer may then be written 
cuz dU? (1) 
ex dx eu? 
the boundary conditions at the moving boundary &% = 0 (say) and at the 
outer edge of the boundary layer being, respectively, 
u(x, 0) 4, lim u(x, ys) U(x). (2) 
Since the velocity on cach closed streamline returns cyclicly to its initial 
value, we have the additional condition that, for each closed streamline, 
u(O, ys) = u(27, (3) 
‘The assumptions that underlie the series expansions of the velocity 
in the boundary layer are best introduced by referring to the whole motion 
of which the boundary layer is part. A parameter y is presumed to be 
formed from the lengths and velocities that enter into the data for the 
complete motion. For a certain value of y (y = 0, say), the motion is 
supposed such that there is no shear layer at the boundary considered, the 
velocity of the inviscid motion exactly matching that of the boundary. 
‘The motion is then considered for the small, non-zero values of y for which 
the velocity of the inviscid motion at the boundary is slightly different from 
that of the boundary. For these values of y it is supposed that LU’? may be 
represented in the form 
UXx) = 1+ YO,(x)y". (4) 
nod 
A similar expansion for u? may then be expected to hold in the boundary 
layer. Accordingly we write 


w(x, = 1+ q,(x,d)y". (5) 
n=l 
Since the boundary laver equation involves wu linearly we also write 
12 
= 41+ Saxby") 1+ Pay", (6) 
nol n=] 


where the functions P,, are polynomials in q;, ‘The viscous term 
of the boundary layer equation is then 


=a t+ ~ +R, )y" 7 

n—-1 q, 

where R,= D> as - (8) 
r=1 ous 


The remainder terms R, depend only on the lower order coefficients 
1s Yoy-++5 J, Of the expansion for u? and may readily be calculated from (6) 
and (8). In particular, 


) 
Ry = 3% 
(9) 

| 

Ry = Ag?) | 
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The boundary layer equation (1) and boundary conditions (2) may now 
be replaced by a set of equations and boundary conditions for the coefficients 


g,- Thus we have 
| 
dx 
with the boundary conditions, for all n, 
g.(x,0)=0, lim g,(x,4)=0,(x), = 9,27, (11) 


(10) 


yon 
Since the velocity gradient udéu/cy must vanish at the outer edge of the 
boundary layer, we also have that 


lim im (x, p) = 0. 


At first sight, it appears that these equations enable all the coefficients q,, 
to be calculated in turn, thereby determining the motion in the boundary 
layer. In fact, however, unless the distribution of the external velocity 
round the boundary layer is suitably restricted, none of these equations 
possesses a solution which satisfies its boundary conditions. For, on 
integrating (10) round a closed streamline, we find that 


aE |, q, dx = 0, 

Rd 

xp}, (n > 1) 


whence, on integrating twice with respect to y and using the boundary 
conditions (11) for q,, 


| O, dx = 0, (12) 

0 

|, O,, dx = | R(x, (n > 1), (13) 
From equation (12) we see immediately that the average value of the leading 
term Q, of the expansion for U? must be zero. As regards equation (13), 
each function R,, is defined in terms of the lower order coefficients q,, gg,..., 
g,—-1 and hence is ultimately determined, on solving (10) successively for 
91» In—1» DY the lower order coefficients Q,, Q,,...,Q,,;. Equation (13) 
therefore implies a relation between each of the coefficients O,, (n > 1) of 
the expansion for U? and the earlier coefficients Q,, Q,,..., O,_,. Unless the 
distribution of the external velocity round the outer edge of the boundary 
layer is such that the conditions represented by equations (12) and (13) are 
satisfied, the boundary layer motion with closed streamlines is not possible. 
The implication of these conditions for the inviscid flow will be discussed 
later. For the moment, it is assumed that the inviscid flow is such that 


they are satisfied. 
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In order to solve the equations for the coefficients g,, of the expansion 
for u®, the functions that occur in these equations are expressed in Fourier 
series. ‘Thus, we write 


(x) > O,, 5 ( 1 4) 


where, from the conditions (12) and (13) just derived, 
Op = 0, 


~ 


Ow=| | Roll’) (n>), 
The equations and boundary conditions for the coefficients g, then reduce 


to the equations 


= tmO 


> 


| 
| (15) 
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with the boundary conditions for all n, 


Gun) = 0, lim (16) 


Your 


These equations define q,,,,, in terms of O,,,, and R,,,,.. On solving them and 
substituting in (14) we get 


L r 
) 24/(im) f ( ) 
where 


amy 

Tam =e { diy + 


- - i 
+ eo vin R,,e m) — I, wine ay 
/0 


As may be recalled, O,,,, can be calculated directly from the external speed 
U(x) (see (4) and (14)). Also R,,,,(y%) can be calculated in terms of the 
coefficients J, 1 of lower order in (see (8) and (14)). Equations (17), 
therefore, may be used to calculate in turn the successive coefficients q,, of 
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the expansion for u?. 
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When the coefficients g, have been calculated in this way, the velocity 
component u may be determined. Also, the conditions (13) may be expressed 
directly in terms of the coefficients O,,,,.. ‘These conditions are now seen 
to be sufficient as well as necessary to allow g, to be calculated. Subject 
to all the series used being convergent, the above calculation thus leads to a 
formal solution for the motion in the boundary layer and enables the necessary 
and sufficient conditions for the existence of the solution to be formulated 
in terms of the external velocity. 

It remains to see how these results fit in with the expected behaviour of 
the inviscid motion. When an inviscid motion is bounded by surfaces 
whose position is known, there is in general one disposable constant of the 
inviscid motion for each closed surface. ‘The conditions on OQ,,,, are clearly 
equivalent to a condition for such a disposable constant. ‘hus on the one 
hand the inviscid motion can generally be chosen so that the boundary layer 
motion of the type considered is possible. On the other hand, the existence 
of the closed boundary layer being assumed, the conditions on O.,,,, provide 
just that extra information needed to determine the inviscid motion. 

A formal solution having been obtained for a class of boundary layers 
whose streamlines are closed, it is of some interest to see whether the solution 
shows any features additional to those normally observed for boundary 
layers whose streamlines are open. ‘Iwo such features seem worth mention- 
ing. ‘lhe first is the property that the speed at the outer edge of the boundary 
layer is related to the speed of the boundary. In assuming that the fluid 
just outside the boundary layer moves continuously round the boundary, 
it was already implicit that the circulatory movement of the boundary could 
impart circulation to the inviscid motion. Some such relation was, therefore, 
to be expected. ‘The second feature concerns the velocity profile. Instead 
of tending monotonically to its limiting value in the mainstream, the 
magnitude of the velocity may oscillate infinitely often about this value. 
The underlying reason for this oscillatory behaviour is clear. If the velocity 
gradient ucu cys tended to zero monotonically, then on streamlines 
sufficiently far from the wall viscous action would always tend to accentuate 
the ditference between the speed in the boundary laver and the speed just 
outside it. As a consequence, the speed on such streamlines would not 
return cyclicly to its initial value. Mathematically the oscillatory behaviour 
results from the periodic dependence of uw on arc-distance, and bears an 
analogy to the oscillatory variation of amplitude with depth in the well- 
known ‘skin effect’ for a rapidly alternating electric current. 

It may be noted in passing that the perturbation method used above 
extends to closed boundary layers in which the fluid is compressible. In 
this case the working assumption is that the motion in the boundary layer 
region differs only slightly from a motion with uniform speed and uniform 
temperature. In order for the closed motion to be possible, certain 
conditions must now be satisfied by the density and temperature just outside 
the boundary layer as well as by the speed there; the energy equation now 
giving rise to a set of conditions analogous to those obtained above from the 
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momentum equation. Both the speed and the temperature may now 
oscillate about their values in the mainstream as the outer edge of the 
boundary layer is approached. 


3. 'THE CONDITION ON THE EXTERNAL VELOCITY DISTRIBUTION 


As has been remarked, the condition on the external velocity distribution 
imposed by the closure of the boundary layer provides information needed 
to determine the inviscid motion. ‘here seems little prospect of deriving 
the exact condition in a simple form. Relatively simple approximate forms 
can, however, be derived, and this we proceed to do. 

The condition on the leading term of the expansion for the external 
velocity is given by (12). Correct to order y this condition is equivalent to 


Ude = 2m, (18) 


- 0 
The circulation round the boundary of the inviscid motion is thus specified 
directly. If the inviscid motion is confined by a single closed boundary 
at which this condition holds then the vorticity of the inviscid motion may 
also be specified directly. For, the streamlines being closed, the vorticity 
w is uniform (Batchelor 1956a). Whence by Stokes’s theorem, 

27 Lity 
where 4 is the area occupied by the fluid. A similar result holds if the 
inviscid motion is confined by two closed boundaries and y is a suitable 


(19) 


expansion parameter for each simultaneously. 

The conditions on the higher order terms of the expansion for the 
external velocity are defined by (13) together with the solution, now known, 
for the motion in the boundary laver. On substituting for R, and R, and 
using the governing equations for qg, and qs, it may be shown that 

“0 Jo Jy /0 
Then, when g, is replaced by the explicit expression obtained for it, this 
relation becomes 


| (U3-U)dx = - ¥ Oy 


= ) 


m 


22(l—m) 


where the overbar denotes a complex conjugate. ‘Ihe condition that the 
boundary layer should close is thus expressed correct to order y? as an 
explicit condition on the external speed. 
Note that to order y?, 


(U-1)dx = 


~ 0 


The circulation imparted to the inviscid flow is thus less than the circulation 


= 0. 


~ 0 
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When the coefficients g,, have been calculated in this way, the velocity 
component u may be determined. Also, the conditions (13) may be expressed 
directly in terms of the coefficients O,,,,.. “Uhese conditions are now seen 
to be sufficient as well as necessary to allow q,, to be calculated. Subject 
to all the series used being convergent, the above calculation thus leads to a 
formal solution for the motion in the boundary layer and enables the necessary 
and sufficient conditions for the existence of the solution to be formulated 
in terms of the external velocity. 

It remains to see how these results fit in with the expected behaviour of 
the inviscid motion. When an inviscid motion is bounded by surfaces 
whose position is known, there is in general one disposable constant of the 
inviscid motion for each closed surface. ‘The conditions on Q.,,,, are clearly 
equivalent to a condition for such a disposable constant. ‘Thus on the one 
hand the inviscid motion can generally be chosen so that the boundary layer 
motion of the type considered is possible. On the other hand, the existence 
of the closed boundary layer being assumed, the conditions on Q,,,,, provide 
just that extra information needed to determine the inviscid motion. 

A formal solution having been obtained for a class of boundary layers 
whose streamlines are closed, it is of some interest to see whether the solution 
shows any features additional to those normally observed for boundary 
layers whose streamlines are open. ‘I'wo such features seem worth mention- 
ing. ‘The first is the property that the speed at the outer edge of the boundary 
layer is related to the speed of the boundary. In assuming that the fluid 
just outside the boundary layer moves continuously round the boundary, 
it was already implicit that the circulatory movement of the boundary could 
impart circulation to the inviscid motion. Some such relation was, therefore, 
to be expected. ‘lhe second feature concerns the velocity profile. Instead 
of tending monotonically to its limiting value in the mainstream, the 
magnitude of the velocity may oscillate infinitely often about this value. 
The underlying reason for this oscillatory behaviour is clear. If the velocity 
gradient ucu cys tended to zero monotonically, then on streamlines 
sufficiently far from the wall viscous action would always tend to accentuate 
the difference between the speed in the boundary layer and the speed just 
outside it. As a consequence, the speed on such streamlines would not 
return cyclicly to its initial value. Mathematically the oscillatory behaviour 
results from the periodic dependence of u on arc-distance, and bears an 
analogy to the oscillatory variation of amplitude with depth in the well- 
known ‘skin effect’ for a rapidly alternating electric current. 

It may be noted in passing that the perturbation method used above 
extends to closed boundary layers in which the fluid is compressible. In 
this case the working assumption is that the motion in the boundary layer 
region differs only slightly from a motion with uniform speed and uniform 
temperature. In order for the closed motion to be possible, certain 
conditions must now be satisfied by the density and temperature just outside 
the boundary layer as weil as by the speed there; the energy equation now 
giving rise to a set of conditions analogous to those obtained above from the 
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momentum equation. Both the speed and the temperature may now 
oscillate about their values in the mainstream as the outer edge of the 
boundary layer is approached. 


3. 'THE CONDITION ON THE EXTERNAL VELOCITY DISTRIBUTION 


As has been remarked, the condition on the external velocity distribution 
imposed by the closure of the boundary layer provides information needed 
to determine the inviscid motion. ‘There seems little prospect of deriving 
the exact condition in a simple form. Relatively simple approximate forms 
can, however, be derived, and this we proceed to do. 

The condition on the leading term of the expansion for the external 
velocity is given by (12). Correct to order y this condition is equivalent to 


7 
| Udx=2n. (18) 
0 
The circulation round the boundary of the inviscid motion is thus specified 
directly. If the inviscid motion is confined by a single closed boundary 
at which this condition holds then the vorticity of the inviscid motion may 
also be specified directly. For, the streamlines being closed, the vorticity 
w is uniform (Batchelor 1956a). Whence by Stokes’s theorem, 


G@ = 


2a Lu 
(19) 
where 4 is the area occupied by the fluid. A similar result holds if the 
inviscid motion is confined by two closed boundaries and y is a suitable 
expansion parameter for each simultaneously. 

The conditions on the higher order terms of the expansion for the 
external velocity are defined by (13) together with the solution, now known, 
for the motion in the boundary layer. On substituting for R, and R, and 
using the governing equations for g, and q,, it may be shown that 

| (l 3_ ‘\dx = ly | | | — g(x, ds O(y'). (20) 

Jo /0 
Then, when q, is replaced by the explicit expression obtained for it, this 
relation becomes 


| (U3-U)dx = - 


/@0 Lm +1ym)? ~ 


 l+m = ) 


m=1 
where the overbar denotes a complex conjugate. ‘The condition that the 
boundary layer should close is thus expressed correct to order y* as an 
explicit condition on the external speed. 
Note that to order y?, 
| 
~ 0 


The circulation imparted to the inviscid flow is thus less than the circulation 
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of the motion at the boundary. ‘To order y the circulation imparted is 
the same as that at the boundary, as may be seen directly from (18). 

In several closed boundary layers to which these results apply, there 
is only one component in the Fourier series for O,(x). Condition (21) 


then reduces to 
2a 


dx = (22) 
0 


4, EXAMPLES 

To illustrate the way in which the additional boundary conditions 
just obtained may be used to determine the inviscid motion, two examples 
are considered. 

The first concerns the motion of fluid enclosed within a fixed, nearly 
circular elliptic cylinder which is lined by a moving band. ‘The band 
moves transversely to the generators and is assumed to maintain the fluid 
in a steady rotary motion in which the action of viscous forces is negligible 
save in a closed boundary layer on the band. 

If the elliptic boundary were exactly circular, the fluid would rotate 
as a rigid body and at the boundary we would have u = U=1. In the 
case of a nearly circular elliptic boundary, therefore, the velocity in the 
boundary layer may reasonably be expanded in powers of the small 
eccentricity ein the same way as in the general case the velocity was expanded 
in powers of y. 

Since the streamlines of the inviscid motion are closed, the vorticity 
of the inviscid motion is uniform (Batchelor 1956a). For a given value w 
of this vorticity, the inviscid motion is determined by the condition of zero 
mass flux normal to the cylinder. It is supposed that the cylinder section 
has semi-axes a, ) and is defined in Cartesian coordinates €, 7 by the equation 


The stream function of the inviscid motion may then be written 


= — (= + + constant 23 
whence the speed of the inviscid motion at the boundary is 
warb? n*\12 
(5 +B) (24) 


To determine the magnitude of the vorticity we use the further boundary 
condition imposed by the requirement that the boundary layer should 
close. As a preliminary we note that to order e, 

U? = 14+(O,)+ } cos 2x)e. 
Moreover, from the condition (12) on the leading term O,(x), Q,) must 
be zero. ‘The Fourier expansion of O,(x) in this example thus contains 
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only one term. ‘The appropriate condition for the velocity at the outer 
edge of the boundary layer is, therefore, that given by (22). On sub- 
stituting for U from (24) and integrating round the boundary, it follows 
immediately that correct to order e%, 

(a? + b?)8 
a®b?{2(a® + + (a2 — 


Whence, after some reduction 


= 


_ 

(ab) 
This relation determines the vorticity and thence the stream function, of 
the inviscid motion correct to order e*. ‘To this order, the vorticity is 
apparently inversely proportional to the square root of the area occupied 
by the fluid. 

Our second example concerns the steady streaming of unbounded fluid 
past a rotating cylinder. A series of photographs which illustrate this 
motion for various values of the peripheral speed u, relative to the mainstream 
speed U,, has been published by Prandtl & Tietjens (1936, plate 7). When 
uy U, is small the boundary layer at the cylinder separates and a wake 
forms. As u)/U, increases, the size of the wake diminishes and ultimately, 
for uo/ U’,, greater than 4, the boundary layer appears to be wholly attached 
to the cylinder. We assume, then, that for sufficiently large values of uy/U, 
the boundary layer remains attached to the cylinder. 

The inviscid motion is then known apart from the value of its circulation 
I about the cylinder. At the cylinder, the speed of the inviscid motion is 
given by 


{1+ O(e*)}. (25) 


(67) 


> 


a, U = +2U, sinx, (26) 


where a denotes the cylinder radius and « the inclination of the outward 
normal to the cylinder to the velocity of the mainstream. 

If the fluid were at rest at infinite distance from the cylinder (U, = 0), 
there would be no boundary layer at the cylinder and I’ would be equal 
to 27a. For values of U, small compared to uy, we may therefore expect 
that U* can be expanded in powers of U,,/u) in the same way as it was 
formerly expanded in powers of y. 

‘lo determine I’ we again use the additional boundary condition imposed 
by the closure of the boundary layer at the cylinder. As in the previous 
example, the appropriate condition is that given by (22). On substituting 
for U from (26) we get 

4 \) 
= 2naw, 11-3 +0() (27) 


4 
ty 


} 
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, L being the lifting 


The lift coethcient for the cylinder, = 
force and p the fluid density, is therefore given by 


C.= = + (28) 
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5. ‘THE CLOSED BOUNDARY LAYER WITH CONSTANT PRESSURE 

So far the speed of motion round the outer edge of the boundary layer 
has been supposed to be non-uniform, the speed of motion at the boundary 
being uniform. It is now supposed that the speed u, U of motion round 
the outer edge of the closed boundary layer is uniform, the speed u f(x) 
of motion at the boundary being non-uniform. Boundary layers of this 
kind are of interest because of the simplicity of the condition they impose 
on the external flow. 

Since the pressure gradient of the external flow vanishes, the equation 
of motion in the boundary layer may be written 


(29) 
ox us? 
On integrating round a closed streamline we get that 
q2 
dx = 0. (30) 
0 
Whence, on further integration, 
1 
C2 = — ax. 31 
(31) 


This is the condition that must be satisfied by the external speed if the 
boundary layer is closed and the pressure gradient is zero. Stated simply, 
the condition is that the speed of the external How must be equal to the 
root-mean-square speed of the closed bounding surface. 

It is rather interesting that this result is unimpaired if the boundary 
layer equation is replaced by the approximate equation 

Cu 

Ox oy? 
Inasmuch as lis a (root-mean-square) average round each streamline of 
the quantity « which it replaces, this equation may well yield a solution 
more accurate than would normally be expected from a iinearized form of 
the boundary layer equations. 

As an example in which condition (31) may be used, we may mention 
the steady two-dimensional motion of fluid inside an infinite circular 
cvlinder which is rotating about its axis, a fixed sheath shielding a segment 
ot the boundary. ‘The inviscid motion is then a rigid-body rotation. ‘Thus 
the boundary layer on the sheath and exposed cylinder surface is closed 
and subject to zero pressure gradient. Hence from (31) the angular velocity 
of the rigid-body rotation Q is related to the angular velocity of the cylinder 
Q, by 


(32) 


Q= (33) 


where x 1s the angle subtended at the axis by the unshielded segment of the 
evlinder surface. 

The author first obtained the result (31) after seeing an early draft of a 
paper on closed flows by Dr Batchelor, and it has been reported in the 
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published version of that paper (Batchelor 1956a). This relation and a 
i weaker form of (21) have been independently obtained by Feynman & 


: - Lagerstrom, who reported their conclusions in a paper presented at the 
; @ 9th International Congress for Applied Mechanics, Brussels, in September, 
> 1956. 

_ ‘This example is also the one mentioned in the Introduction as having 
y been considered by Squire. Squire used a Rayleigh-type approximation 


in which the boundary layer equation was, in effect, replaced by 


A a9 
Ou 
U — 
ox oy” 
) | v being distance normal to the boundary. ‘This equation yields the linear 
relation 
= a. (34) 
) 
which is, of course, correct when x = 0 or 27 but elsewhere gives a value 
for 2 which is too small. ‘The greatest error in 2 occurs when « = 37, 
) the result then being too small by a factor }. 
The author is very grateful to Dr I. Proudman for his encouragement 
: and advice and to Dr G. K. Batchelor whose interest in motions with closed 
° streamlines gave the initial stimulus to this work. 
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A general solution of the equations of hydrodynamics 


By L. M. MILNE-THOMSON 
Brown University, Providence, Rhode Island 


(Received 9 October 1956) 


Consider a fluid of constant viscosity ». If the fluid is incompressible, 
We suppose it to move in a conservative field in which the force per unit mass 
is —~VQ. If the fluid is compressible, we suppose the field of force to be 
absent. Both cases are included if we write V = epQ, where « = 1 for an 
incompressible fluid, and « = 0 for a compressible fluid. ‘The equation of 
continuity 


+ V. (pq) a (1) 
can be satistied identically by 
p=N*x, pq = —V(ex/et), (2) 


where y is an arbitrary (differentiable) scalar function of position. 
The equation of motion is (Milne-Thomson 1955, §§ 19.02, 19.03) 
p dq/dt = —-VV+V.9, (3) 
where the stress tensor is 
—(p+ (+) 
and / is the idemfactor or unit tensor of the second rank. 
Now (3) can be written 
pdq/dt V1). (5) 
Using (1) we find that (Milne-Thomson 1955, § 2.34) 
p dq dt = (pq)/ct + V.(pq;q) = V.(— + pq; q) 
and therefore (5) becomes 
Dp af I(6*x ot? 1’) - pq: ==) 
which can be satisfied identically by 
= (6) 
wherein ‘Vis an arbitrary tensor of the second rank (Milne-Thomson 1942; 
that V.(Va‘4'a¥) =0 is a consequence of Theorem I of this paper). 
Equations (2) and (6) furnish the density, velocity, and stress distribution 
in terms of an arbitrary function and an arbitrary tensor of the second rank. 
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A boundary layer theorem, with applications to 
rotating cylinders 


By M. B. GLAUERT 
Department of Mathematics, University of Manchester 


(Received 9 November 1956) 


SUMMARY 

If, in a given solution of the boundary layer equations, the 
position of the wall is varied, then additional solutions of the 
boundary layer equations may be deduced. ‘The theorem 
considers the nature of such solutions, for the general case of 
time-dependent three-dimensional compressible flow. 

Applications of the theorem. arise in several different fields, 
and it is shown that useful quantitative results can often be 
obtained with the minimum of calculation. In this paper, chief 
attention is focused on the case of a rotating circular cylinder, 
and explicit formulae are developed for the skin friction, valid 
for sufficiently low rotational speeds. ‘The important results 
which the theorem gives for slip flow have been noted by previous 
authors, and only a brief discussion is given here, but certain 
extensions to these previous treatments are made. Other 
applications of the theorem are briefly mentioned. 


1. INTRODUCTION 


‘The theorem concerns the additional solutions of the boundary layer 
equations that may be deduced from a previously known solution, by 
considering the body surface to be at z= {(x,y,¢) instead of at z= 0, 
where (x,y) are orthogonal curvilinear coordinates on the body surface 
and ¢ is the time. For the case of steady two-dimensional incompressible 
flow, the theorem was given by Prandtl (1938), though he does not appear 
to have made any applications of the result. More recently Nonweiler 
(1952) effectively uses a limited form of the theorem for steady two- 
dimensional compressible flow, in a paper concerned with the effect of 
slip in a laminar boundary layer. He introduces the concept of a ‘ plane 
of zero-slip’, at a small distance below the actual surface, and applies the 
usual no-slip boundary conditions there. He deduces results equivalent 
to those following from the theorem, though his treatment is somewhat 
imprecise. Mangler (1956) has attempted to extend Nonweiler’s work 
to steady three-dimensional flow by employing a three-dimensional form 
of the theorem, but the boundary layer equations he uses are inadequate, 
as all curvature terms are omitted. 
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In this paper the theorem is proved for the general case of time-dependent 
three-dimensional compressible flow. ‘The form of the new surface 
z= {(x,y,t) is shown to be arbitrary, with the proviso that ¢ shall not 
vary so rapidly that extra curvature terms have to be included in the 
boundary layer equations, and the conditions for this to hold are investigated. 
It might be anticipated that difficulties would arise in attempts to apply 
the theorem when ¢€ 0, since in this case the new solution is not defined 
for ¢< zs <0. However, a given original solution of the boundary layer 
equations can be extended for a certain distance across the surface 3 = () 
by analytic continuation, and it will be found that in practical applications 
¢ is only a small fraction of the boundary layer thickness. Numerical 
results will be obtained in terms of values and derivatives of the dependent 
variables at s = 0 in the original solution, and thus the theorem will be of 
use for both positive and negative values of €. 

Fields of application of the theorem include flow past a rotating cylinder, 
slip tlow, flow over a body covered in whole or part by a moving belt or a 
liquid film, and flow over a body of liquid. Some consideration is given 
to all these cases in this paper. 


2. ‘THE ‘THEOREM 
The boundary layer equations for unsteady compressible three- 
dimensional laminar flow can be written as follows, in terms of orthogonal 
curvilinear coordinates (x,y, 3), where z = 0 is the body surface: 


p Dt ex cv oz 
Du puv ch, pi? chy le cu\ | 
Dz pu* ch, puv ch, 1 


Dp Cu 2 cu 2) , 


Here A, dx and h,dy are elements of !ength in the coordinate directions, 
h, and h, being tunctions of v and y, (#, 7, «w) are the velocity components, 
wand & are the coefficients of viscosity and thermal conductivity, and ¢ is 
the time. ‘The pressure p, the density p and the temperature T are also 
connected by the equation of state. The pressure does not vary across 
the boundary layer, and is known from conditions in the external flow. 

Suppose one solution of these equations, Solution 1, is available, where 
the variables take values (u),7,,%@,,p,, 7) which are known functions of 
(v,¥,3,¢). They will satisfy certain conditions at the surface s = 0, and 
at the outer edge of the boundary layer, which may be considered to be 
at s = « as far as the boundary layer equations are concerned, all except 


| 
a 
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w, wil! take prescribed values. Consider now a second set of expressions 
(3, Py, with the properties 

u(x, y, 2, t) = u,(x, y, 2+ ¢, 2), 

U(x, 2, t) = 0,(x, y, 2+ ¢,t), 


pAx, t) = p,(x, y, 2+ ¢, 2), 


(2) 
y, = T(x, y,2+¢,2) | 

= feo, - 
WAX, V, = {WwW = | 

ct h, Ox h, ey) tar J 


where ¢ = {(x,,¢). Our theorem consists of the assertion that the 
quantities (2) also satisfy the equations (1), and so give a second solution, 
Solution 2, of the boundary layer equations. At the outer edge of the 
boundary layer, Solution 2 will satisfy the same conditions as Solution 1. 
The value of w will be different, but this is not a quantity which can be 
prescribed in a boundary layer problem. At the surface, s = 0, Solution 2 
will satisfy quite different conditions from Solution 1, and it remains to 
be investigated whether they will be appropriate in any practical application. 

‘The general form of the expressions (2) is in accord with the concept, 
mentioned in § 1, of taking a new surface at s = Cin Solution 1. The extra 
contribution —c¢/ct to w, results from the relative motion of the old and 
new boundaries, and the contribution —(u,A,)e¢ ex ey arises 
since u, and 7, have components normal to the new surface. In boundary 
layer flow, w is small compared with u and z, so there are no corresponding 
additions to wu, and 7. 

However, the truth of the theorem is most concisely and conclusively 
demonstrated by direct substitution in the boundary layer equations (1). 
From equations (2), 


(3) 
os os 

Dt cl h, Ox ly 


(4) 


and similarly for p,, and 7T,. Also, 


uy) t 2 Us) 


: 
4 
1 
ai 
rq 
= 
) 
Cu, cu, u,feu, cu, v, (Cu, ou, eC 
ct ds ot h,\ ex OX hy oy oz 
et) ox hy dy] az 
ot 6h, ox Oy dz 
D 
l i's 1 uy) 
hyh, ex ey oz hyhy| ox Oz ox 
ov loz \ds dz dx hy dz oy} 
t Ox hy oz 
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It is thus clear that equations (1) and the equation of state are satisfied, 
since p, h, and A, are not affected, and so the theorem is proved. 

We should expect there to be some limitations on ¢ for the theorem 
to hold. In boundary layer flow, ¢?u/é3? is O(6-?), and no other second 
derivative of u exceeds O(é-!), where 6 is the thickness of the boundary 
layer, and the dimensions of the body are considered to be O(1). This 
must continue to apply, or extra viscous terms will have to be included 
in equations (1). Now 0?u,/éx? has terms 


and 


and hence we must have 


> 


r 


| 


= = = O(1). (6) 


Consideration of other derivatives leads to similar conclusions. Clearly 
€ = O(8), or Solution 2 will lie entirely outside the original boundary 
layer. We see, therefore, that the typical distance for a change in ¢ must 
not be less than O(6!?), i.e. a distance intermediate between the length of 
the body and the thickness of the boundary layer. 

The direct physical interpretation of the theorem, as given by 
equations (2), is that a given solution may be displaced through a 
distance —C€ in the z-direction, without affecting its validity. This is 
not exactly the same as our previous idea of taking the body surface to 
be at x = ¢ instead of at z = 0, as in the latter case a change of coordinate 
system is involved, the new coordinates being along and perpendicular to 
the new boundary. However, it would seem that on the boundary layer 
approximation the two physical processes are indistinguishable. If the 
change of coordinate system is not to affect the values given in equations (2), 
and if the curvature of the boundary is to be such that the solution i 
applicable to bodies of small or zero surface curvature, it can readily be 
shown that ¢ must satisfy the conditions of equation (6), precisely as before. 
‘This gives further confirmation of the equivalence of the two physical 
interpretations of the theorem. ’ 


3. ROTATING CYLINDER 
Practical applications of the theorem can be found when ¢ is a small 
fraction of the boundary layer thickness 5. For the simple case of steady 
two-dimensional incompressible flow, we may take #, = 1, and equations (2) 
reduce to 


\ 


u(x, 2) = u,(x, 2+), w(x, 3) = w,(x, 3+ ¢)—Uu,(x, 5+ C) » 


where ¢ = (xv). In terms of the stream function &, equations (7) are 
equivalent to 


(x, = (x, (8) 


\ 
| 
: 
| 
} 
a" 


tv 


A boundary laye: theorem, with applications to rotating cylinders 93 


in which form the theorem was noted by Prandtl (1938). By use of Taylor’s 
theorem, we now obtain 


OY; 7} 
be = f+ +.. 
Ou | 
(9) 
ty 
Ole Ou, , | 
+.... | 
CF C2 oz? J 


At s = 0, Solution 1 will normally satisfy Y, =u, = 0. Also 
T 
(=) = +, (10) 
where 7 1s the skin-friction, and 


_ 


oz 


from the boundary layer equations, where U(x) is the velocity outside the 
boundary layer. ‘hus if (¢/6)? is negligible compared with unity, the 
conditions satisfied by Solution 2 at s = 0 are 

= 0, tte = Cr,/p, T, = T,—CpU(dU/dx). (12) 


Consider the flow past a circular cylinder, rotating so that the surface 
has a small tangential velocity a. If Solution 1 gives the flow with zero 
velocity at the cylinder surface, then by (12) we obtain the required 
Solution 2 by choosing 


= (13) 
(14) 


‘Thus if 7, has been determined, perhaps by means of a Karman—Pohlhausen 
calculation, 7, is given at once. It should be noted that the external flow 
U(x), used in (14) and to obtain 7,, must be that round the cylinder when 
it is rotating, and the determination of this flow is not a purely boundary 
layer problem. ‘The condition for (12) to hold, that ¢/6 is small, implies 
that a/U must be small. ‘This is clear physically, since the new surface 
in the Solution 1 flow must be at the point at which u, = a. As discussed 
in §1, equation (14) applies whether a is positive or negative. 


Separation and stagnation points : 


According to equation (13), ¢ becomes large when 7, becomes small, 
and so our solution ceases to be valid, whatever the value of a. This will 
occur in the region of the front stagnation point and at separation. The 
question of what constitutes separation on a rotating cylinder is one of 
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some complexity, which we shall not embark on here. We note that, 
according to (14), 7, becomes infinite when 7, = 0. However, if we retain 
the terms in @ in (9), the boundary condition on w Is satisfied if 


ap. = \@pU(dU/ dx), (15) 
and with ¢ chosen to satisfy (15) we obtain 
t, = {ri (16) 


‘This expression for 7, is equivalent to (14) when 7, is not small, but does 
not become large near separation. If a@ is small, ¢ is small for all values 
of 7,, and the value of ys, at s = 0 never exceeds O(a U)??. ‘Thus (16) is 
applicable even near the separation point in Solution 1. 

Consider the cylinder to be rotating so that the upper surface is moving 
in the same direction as the stream. ‘Then a is positive on the upper surface 
and, according to equation (16), 7, has a finite value when 7, = 0 (since 
dU’ dx is negative here). ‘Vhe validity of Solution 1 ceases at this point, 
and so no further deduction of Solution 2 can be made, but presumably 
t, falls steadily to zero as xv increases further. On the lower surface a is 
negative, if xv is still measured in the flow direction, and 7, = 0 when 
7 = 2aupl (dl dx). ‘Vhe solution cannot be extended further by use of 
the theorem, as there is no point in Solution 1, or its continuation below 
x = 0, which could be taken as the new position of the surface. In any 
case, it is natural to consider 7, to be zero beyond this point. 

Near the front stagnation point, where l’-— a, there is no point in the 
boundary layer at which uw, = a, and so the theorem must cease to apply. 
However, the value of 7, given by equation (14) has a finite limit at the 
stagnation point, and we can show that this value is correct. Consider 
the linear approximation to (9), 

= Wy + Cuy, Uy = U, + C(eu,/ez). (17) 
If (eu, es). = a, the boundary conditions at z = 0 are satisfied exactly, 
and the boundary layer equations are satished if terms quadratic in ¢ are 


“ \ 
eu,\ | feu, 

= 4,4 af = (2) (18) 
\ 7 2=0 


In this form it is clear that, in the boundary layer equations, the terms in a? 
remain finite as v—-0, and will be negligible if @ is small, and so the 
equations are adequately satisfied. Equation (14) follows as before, and 
so is valid even in the neighbourhood of the stagnation point. Incidentally, 
Rott (1956) and Glauert (1956) have shown that, for a linear external 
velocity distribution U’ = cx, equation (18) gives the exact flow for a surface 
velocity a of any magnitude. 


neglected. “Thus 


Torque on cylinder 

We now attempt to make a quantitative estimate of the torque G ona 
slowly rotating cylinder. ‘The first thing to be decided is the form of the 
velocity distribution U(x) at the edge of the boundary layer. Experiments 
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by Reid and by Betz quoted by Goldstein (1938, p. 546) show that no 
, q significant lift is observed if a — }l, where Uy is the velocity of the 
oncoming stream, and so we shall assume that U(x) remains the same as 
when the cylinder is not rotating. In accordance with the values given 
by Hiemenz (Goldstein 1938, p. 150), we choose 


U/U, = 1:81380 —0-270568 — 0-047086, (19) 

where 4 is the angular coordinate round the cylinder. For a non-rotating 
" cylinder, the skin-friction is now given by the series due to Howarth (1935), 
; : with additional coefficients calculated by Ulrich (1949), as 
¥ 

RI? 42580 — —0-28765 + 0-01906" — 0-00046", (20) 

} Puy 
° ; where the Reynolds number R = U, dv, d being the cylinder diameter 

] and v the kinematic viscosity. 

On the basis of (19) and (20), the value of (7, —7,)/a as given by equation 
. : (14) can be tabulated as « function of 4, and the results are shown in figure 1. 
: : Over the front part of the cylinder, the skin-friction produces a decelerating 
~ a torque, tending to reduce the rotation, but beyond the velocity maximum 
f there is an accelerating contribution to G/a, which according to (14) is 
» | logarithmically infinite. ‘There is thus a tendency to autorotation. 
ab 

6 x| 
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Figure 1. Increment of skia-friction (7.— 7,)0n a circular cylinder rotating with surface 
velocity a in a stream Up, at an angle @ from the front of the cylinder. 


te 


very small (equation (14)), ———- -——-—- a/U, = 0-025, 
—— —— a'U, = 0°1, x limit of validity. 
] However, this result is valid only for vanishingly small a; for larger 
‘ values of a, equation (16) must be used in place of equation (14), in the 


separation region. ‘The values given by (16) cannot be represented by a 

single curve for all values of a. In figure 2 the values of 7, are shown 

directly for a/U, = 0-025 and a/U, = 0-1. Here equation (14) has been 
a used to obtain the values up to the velocity maximum (where 7, = 7,), and 

equation (16) from there on. The corresponding values of (r,—7,)/a@ are 
$ also shown in figure 1. 


ak 
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We now calculate the torque by integrating the skin-friction over the 
cylinder, as far as the zero of 7, on the lower surface and the zero of 7, on 
the upper. As discussed above, these are the limits of validity of our solution. 
We obtain the following results. 


a = G 

When - = 0-025, Ri? = — 0-0088. | 
(21) 

ad G 

When 0-4, R pled? 0475. 


‘The negative sign indicates that the torque G acts to decrease the rotation. 
Certain sources of error in the values of equation (21) may be noted. On 
the upper surface, 7, must fall to zero over a certain distance beyond the 
point at which 7, = 0, and this will give an additional accelerating torque. 
However, it would appear from figure 2 that the final values of 7, are probably 
overestimates, due to the increasing inaccuracy of our approximate solution, 
and also there may be a significant torque over the separated region at the 
rear of the cylinder. ‘These two effects will produce an additional 
deceleration. Further errors will have arisen from our choice of a 
symmetrical velocity distribution, and in sum there seems little justification 
for modifying the values of equation (21), though no great confidence can 
be placed in their accuracy. 


Figure 2. Skin friction 7, on a circular cylinder rotating with surface velocity a in 
a stream Us, at an angle @ from the front of the cylinder. 
——— —— a/U, = 0 (no rotation); —— - —— - a/U,) = 0-025 ; 
—— aU, = 01. 


Finally, a rough calculation shows that if a Us is less than 0-001, G is 
positive and so tends to produce autorotation. On the basis of this estimate, 
a cylinder falling freely with velocity Uy, should rotate, with a surface 
velocity 0-001 

For a cylinder making rotational oscillations with frequency 2, our 
solution will continue to apply provided ¢c{/ct is not too large, so that, 
from equation (2), the boundary condition on w, is adequately satisfied. 
It is easily shown that this requires that Qd/U, shall be small compared 
with unity. For a cylinder making translational oscillations (Glauert 1956), 
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the torque is expressible as a power series in Qd/U,. No doubt the same 
is true in the present case, the theorem enabling us to estimate the first 
term of the series. 


4. SLIP FLOW 


For the high speed flow of a gas of low density, the boundary conditions 
at the surface may be modified to permit a slip velocity and a temperature 


jump. In the notation of §2, the boundary conditions at z = 0 for three- 
dimensional flow are 
u = Lou/oz, 
v = Lov/dz, (22) 
w= 0, 


T-—Ty = L*0T/dz, 

where 7’ is the temperature of the solid surface and L and L* are lengths 
of the order of the molecular mean free path. The basis of these conditions is 
discussed by Nonweiler (1952), who gives an extensive treatment of the 
two-dimensional case, using what is in effect a form of our theorem for small 
¢. We shall show briefly how the theorem may be applied in the general 
three-dimensional case, as was attempted by Mangler (1956), who, however, 
omitted all curvature terms from the boundary layer equations. 

Use of ‘Taylor’s theorem in equation (2) shows that, if (£/5)? is negligible 
and the flow is steady, 


Uy = Uy +Cdu,/dz, ) 

Ug = U, + 

(23) 
Cw, uy vy 

ax hy Oy” J 


Solution 1 will normally satisfy uv, = 7, = w, = 0 at z = 0. ‘The conditions 
(22) on u, v and w will then be satisfied by choosing 
f= (24) 
since, from the first of equations (1), ow,/éz = 0 when z=0. For the 
temperature, it follows from (22), (23) and (24) that, at z = 0, 
Ty = (L- L*) 0T,/oz. (25) 
The following expressions for 7, and 7,,, the components of the skin-friction, 
and q = (kcT/éz),_», the heat-transfer rate at the surface, are now readily 
obtained with the use of Taylor’s theorem and the values of equations (1) 
at s = 0. 


L ep 
= + h, ax 
_ L op (26) 
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Thus, corresponding to a known Solution 1, we can find a Solution 2 
applicable to slip flow, the variation of temperature and heat transfer over 
the surface being given by equations (25) and (27). ‘The extra contributions 
to the skin-friction are given by equations (26), without even requiring a 
knowledge of the Solution 1 values. ‘There are no difficulties either at a 
rounded nose or at separation, since 1/4 remains small. At a sharp leading 
edge the solution ceases to be valid, but so do the boundary layer equations 
themselves. 

For two-dimensional flow, all these results were obtained by Nonweiler, 
although his expression for w, was incorrect. On the assumption that the 
flow is incompressible, Lin & Schaaf (1951) gave the results for the special 
cases of flow over a flat plate and near the stagnation point on a body of 
revolution. 

Difficulties occur in practice. If the wall temperature is to be constant, 
or the rate of heat transfer is to be zero, it is insufficient to employ a Solution 1 
with the same property. As shown by (25) and (27), the required Solution 1 
must satisfy a slightly modified temperature condition, and hence (7,.), and 
(z,,), will differ from their values in the corresponding no-slip case by terms 
proportional to L. Nonweiler makes several attempts to estimate the 
effects of such modifications, 


5. FURTHER APPLICATIONS 

Brief consideration will now be given to certain other classes of problem 
in which the theorem might have useful applications. In all cases the 
velocity distribution outside the boundary layer must be specified, and this 
may be difficult in some instances. 

For the flow past a body covered in whole or part by a moving belt, 
the theorem is applicable if the belt velocity is small compared with the 
stream velocity, the treatment following that of §3. As is clear from 
equation (23) the theorem gives a solution only if the belt moves in the 
same direction as the fluid near the surface, so this application is effectively 
restricted to two-dimensional motions. If the belt extends over only part 
of the body, difficulties arise when considering the flow near its ends, 
as ¢ would have to change more rapidly than permitted by the theorem. 

For the flow past a body covered with a liquid film, Solution 1 would be 
computed assuming the film to be frozen. The skin-friction induces a 
motion in the film, and the effect of this motion on the flow over the body 
could be found from Solution 2. Since the liquid moves in the direction 
of the stress, the theorem will be applicable for a general three-dimensional 
body. 

For the flow past a body of liquid, Solution 1 would be calculated 
assuming the body to be solid. If the internal motion of the liquid could 
then be estimated, its effect on the air boundary layer would be given by 
the theorem. Here the applications will usually be confined to two- 
dimensional or axi-symmetric motions, 
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REVIEWS 


Theory of Combustion Instability in Liquid Propellant Rocket Motors, 
by Lurci Crocco and Sin-I CHENG. London: Butterworths Scientific 
Publications ; New York: Interscience Publishers Inc.; 1956. 
200 pp. 37s. 6d. or $5.25. 

There are some surprises in store for those who have thought of rockets 
as rather superior squibs. It is perhaps no surprise that their history 
starts in thirteenth-century China, where they were used not only as devil- 
scaring fireworks but also as incendiaries. ‘That they were used in the 
eighteenth century by the Indians to break up British cavalry charges 
appears to be substantiated, for it was from India, at the close of the 
eighteenth century, that Sir William Congreve and others brought the 
idea of the rocket as an artillery weapon. Nelson bombarded Copenhagen 
with rockets. Francis Scott Key watched the bombardment of Fort 
McHenry in 1814 from a British ship, and we suppose that rockets were 
used in this action because in the Star Spangled Banner, which he wrote 
the next morning, he refers to ‘‘ the rockets’ red glare’”’. In all of these 
rockets and in the short range artillery rockets used in the Second World 
War, solid propellants (gunpowder, cordite, etc.) were used to generate 
the hot gases. 

The development of rockets using liquid propellants was started in 
1920 by R. H. Goddard in the U.S.A. and by the Society for Space Travel 
in Germany in 1927. Work on liquid-propelled rockets was done at 
Peenemunde in Germany, in California, and at Ministry of Supply Establish- 
ments, during the war. ‘l’oday liquid propellent rockets are being used for 
aircraft propulsion, both as a main power plant and for boost during take-off, 
and for long range missiles, and will be used for firing the U.S. earth 
satellite. Large liquid rocket motors develop thrusts of 100000 1b and 
more. 

The rocket motor for liquid propellants seems at first sight to be 
comparatively simple. ‘he combustion chamber is a tube with an 
injector plate forming one end and a convergent-divergent nozzle at the 
other. ‘lhe chamber and nozzle may be lined with refractory material or, 
more generally, cooled by the propellant. A large number of small holes 
are drilled in the injector plate, and the propellants are forced through 
them, either by pumps or by applying gas pressure to the fuel tanks. ‘The 
pressure in the combustion chamber generally lies between 300 lb per 
sq. in and 600 Ib per sq. in., and the pressure drop through the injector 
plate is of the order of 100 lb per sq. in. 

‘The hot gases may be generated by the exothermal decomposition of a 
single liquid. For example a monopropellant such as hydrogen peroxide 
decomposes with the evolution of heat into hydrogen and steam. Alter- 
natively two reacting liquids, one known as an oxidant and the other as tuel 
may be used. Oxidants commonly used are liquid oxygen and nitric acid 
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Fuels include alcohol, and kerosine. Factors which affect the choice of 
propellants include availability (or price), ease of handling, toxicity, and 
the performance. ‘lhe performance is specified in terms of the specific 
impulse, that is, the thrust produced per unit flow rate of total propellant. 
Specific impulses vary normally from 200 to 250 lb of thrust per Ib per 
second of propellant flow. A rocket developing a thrust of 100000 Ib 
therefore uses fuel at an alarming rate, between 400 and 500 Ib per second. 
At this rate an ordinary motor car petrol tank would empty in about one 
fifth of a second. Large valves and pipes of high capacity are required 
for the propellant feed. ‘lurbo-pumps developing several thousand 
horsepower may also be required. ‘The rate of propellant flow to the 
injectors has to be controlled. Systems have been devised using constant 
injection pressure, as may be obtained with gas pressurization, or constant 
feed rate, as may be obtained with a displacement pump or cavitating 
venturi. Combustion chamber pressure may also be used to control turbo- 
pump speed or feed valve openings. In addition, for bipropellant systems, 
the mixture ratio needs to be controlled. ‘The propellant feed system 
therefore tends to become somewhat complex on the larger motors. 
Reliability is an absolutely essential feature of such systems. In the event 
of loss of flame a large pool of explosive propellant mixture could rapidly 
build up inside the rocket with every possibility of explosion. 

The hazards of starting a rocket motor are considerable. Ignition is 
accomplished in bipropellant systems by an igniter flare, and in mono- 
propellant systems generally with a catalyst. Ignition must be sure, or 
there is again a danger of collection of unburnt propellant which may 
explode. Another hazard of more subtle origin occurs because the build- 
up of chamber pressure in the starting phase introduces transients into the 
propellant feed system and into the combustion process. From these 
transients instabilities may develop which may end in explosions, or holes 
may be burnt in casings and injector plates, thus leading to the collapse of 
the rocket and to fires as the fuel lines are severed. _ Failures develop very 
rapidly, and the detection and prevention of the spreading of a failure is 
difficult because of the short total times involved (no more than 15 to 60 
seconds for liquid rockets). It has been found that oscillations of pres- 
sure in the combustion chamber and propellant lines may be the pre- 
cursor of failure. 

Practically all of what is known about the design and behaviour of 
rocket motors has been found by trial and error. No one designs an 
injector plate with its multitude of small holes, they develop it. ‘The 
high temperature of rockets (now a white rather than a red glare) 
makes the ordinary difficulties of experimentation .in practical com- 
bustion systems even greater. Measurements are limited by these 
difficulties to determinations of thrust, fuel consumption, chamber and 
other pressures, and some temperatures. Detailed exploration of the 
processes of liquid jet impingement, droplet evaporation, mixing and 
combustion within the chamber, have hitherto been practically impossible. 
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Owing to the large density changes from liquid to burnt gas, simulation 
in the absence of combustion is unsatisfactory. ‘lhe assessment of the 
merits of an injector plate depends on a few overall measurements inter- 
preted in the light of the ‘ feeling’ of the investigator. Over the ex- 
tensive period of development which appears necessary, investigators 
grow up, get better jobs, and in this and other ways any small thread of 
continuity in a series of not very precise assessments is lost. ‘lhis process 
has given us injector plates and rocket motors which work, although it 
has not given us any very precise idea of how they work or how to design 
other plates. ‘hose groups who have spent much time and etfort in 
rocket development have acquired know-how, but not enough knowledge 
to make the procurement of a new, large rocket motor anything but an 
expensive and protracted development process. 

It is against this general background of rocket development that one 
must assess the monograph by Crocco and Cheng. ‘The problem of 
instability in rocket motor combustion is difficult to treat analytically, but 
it is the most expensive unknown in rocket development. Instability in 
the combustion seems, at best, to increase the local heat transfer rates, which 
often leads to burning of injector plates or other walls, and at worst it may 
cause oscillations in pressure jarge enough to lead to explosions. In 
current testing practice any detectable vibration is generally the signal for 
a shut-down tor safety. ‘he authors have treated two kinds of instability. 
‘The first is a low frequency *‘ chugging’ in which chamber pressure and 
propellant feed interact to build up oscillations (of 40 to 120 cycles per 
second) in chamber pressure and propellant feed rate. ‘Ihe second, known 
as ‘ screaming ’, is a high frequency oscillation (of 600 to 2000 cycles per 
second) which is related to the fundamental longitudinal acoustical mode 
in the combustion chamber with both ends closed (so that the frequency 
varies inversely with chamber length). ‘lhe treatment is perforce confined 
to the linear instability of small disturbances. ‘lhe basic assumptions in 
such a treatment are all-important, and the assumptions made by the 
authors about the combustion process in particular are worth noticing here. 

{n liquid rockets propellants are injected through the small holes in the 
injector plate as liquid jets. Some of the holes are drilled at an angle in 
order to make neighbouring jets impinge. Fuel and oxidant jets are inter- 
spersed, and the total and relative rates of injection per unit area of face of 
the injector plate may vary with radial distance from the axis of the com- 
bustion chamber. More fuel (relative to oxidant) is often injected near 
the outer wall to reduce gas temperatures at the wall. In the monograph it 
is assumed that the flow is radially unitorm and substantially one- 
dimensional. It might then have been assumed that, as the propellant 
passes axially along the combustion chamber, it changes discontinuously 
trom liquid to vapour and then mixes and reacts to form burnt gas. Such 
an assumption would have made possible an allowance for any delay in 
reaction due to imperfect mixing. Actually the authors prefer to assume 
that the liquid volume is negligible so that the chamber is filled with gas, 
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and that immediately after the liquid evaporates it burns, the time spent 
in the unburnt gas phase being negligible. The time between injection 
and combustion of an elementary mixture is called the time lag. It may 
ditfer with differing propellant elements, and is made up of a portion which 
is insensitive to chamber pressure, this being the time for jet impingement 
and droplet formation, and a portion which is sensitive to variations in 
chamber pressure, being the time for evaporation and combustion. ‘The 
relation between the rate of burning and the rate of injection of a propellant 
element is expressed by the value and rate of change of the combustion 
time lag. 

The authors imagine that the propellants enter the combustion chamber 
as droplets of negligible total volume and, in the steady operating state, 
change to burnt gas. In one case they assume that the change to burnt 
gas of all droplets occurs simultaneously at a specified distance from the 
injector plate, the droplets having passed through a stationary (or re- 
circulating) mass of burnt gas. In this case one may assume that all 
droplets pass along the same path and have the same characteristics, so 
that their time lag is the same. Alternatively there may be a spread in 
time lag due to different droplets having passed to the combustion front 
along different paths, In other cases the ‘ combustion’ is assumed to be 
distributed along the axial length of the combustion chamber. Here 
again different and independent assumptions may be made about the time 
lag distribution. ‘The only question which these assumptions introduce 
(apart from the major difficulties of determining numerical values for time 
lags and combustion distribution) is whether the assumption that all the 
gas is burnt will introduce appreciable error. 

In the treatment of chugging, the frequencies are low enough to permit 
the assumption that the gas pressure is uniform in the combustion chamber 
and varies only with time. Further assumptions are that the temperature is 
constant and that the time lag is uniform, and the authors include an 
examination of the etfect of variations in these assumptions. ‘lhe governing 
equation is then reduced to a mass balance, relating the rate of generation 
of burnt gas to its rate of ejection through the nozzle plus its rate of accu- 
mulation in the combustion chamber. ‘The rate of generation of burnt gas 
is related to the rate of injection of propellants by the time lag. ‘The rate 
of ejection through the propelling nozzle, when upstream conditions are 
oscillatory, is found to depend on the frequency of the oscillations and the 
geometry of the subsonic portion of the nozzle. ‘The solution of the 
resulting equation depends on the dynamics of the propellant feeding 
system. Some relatively simple examples are solved numerically. An 
interesting suggestion which is examined is the use of a servo-stabilizer in 
which a variable capacity of the feeding system is introduced and is suitably 
controlled by the chamber pressure. 

The results of the analysis of chugging should be and have been of 
considerable value. ‘They show the direction in which changes should be 
made to improve stability ; for instance, change of the feeding line length 


1 
: 
1 | 
q 
4 
it 
h 3 
y #4 
n 
d 
i. 
n 
e ‘ 
d | 
n 
Le 
Ag 
nh 
ot 
l- 
ir 
it 
it 
y 
s, 
Rite 
: 


Reviews 


104 


and increase of injector pressure drop both help. Of considerable 
importance is the interaction index m, which occurs in the assumption that 
the combustion rate is proportional to the mth power of the pressure. The 
authors argue that ”, which covers the pressure dependence of the evapo- 
rative as well as of the combustion process, will be of the order unity and 
larger and show that increase of » has a destabilizing effect. 

In the range of frequencies of screaming, the wavelength of standing 
oscillations becomes comparable to the length of the combustion chamber. 
Local pressure variations must therefore be taken into account, so that the 
spatial distribution of combustion and pressure becomes important as well 
as the combustion time lag and the time variation of pressure. ‘The authors 
consider a number of cases, starting with the simplest in which the time lag 
is uniform and combustion is concentrated at one axial position. 

Of major interest is the determination of the boundary conditions at the 
nozzle end of the combustion chamber. In a one-dimensional analysis the 
authors show that the convergent-divergent nozzle behaves more like a 
closed end than an open end to the combustion chamber. ‘The reflections 
are weaker than, and displaced in phase relative to, those at a simple closed 
end. ‘lhe authors show that increasing the length of the subsonic portion of 
the nozzle has a stabilizing etfect. ‘The shape of the converging portion 
also has an effect, the best shape suggested being that which produces a 
velocity distribution linear with the axial length. 

The instability due to transverse (or swashing) oscillations is mentioned 
only briefly. ‘This is unfortunate, since it now appears that this form of 
oscillation is at least as important as the others described. 

One major critical question which may be raised is whether such a 
specialized book is timely. ‘The structure of ideas and the mathematical 
models underlying the analysis rest on too weak a foundation of experiment 
for the volume to be in any sense definitive. Admittedly, one of the objects of 
Agardographs (that is, monographs sponsored by the Advisory Group for 
Aeronautical Research and Development, attached to NA’T'O), of which this 
book is one, is to present current thought and so to stimulate experiment 
and research. On reading this book the rocket researcher might well be 
stimulated to answer experimentally such questions as, what is the value 
of the interaction index, what are the neutral stability points, and what 
is the effect of injector plate design on these and other factors? But could 
not the authors’ ideas have been communicated better in the customary 
medium of papers in journals, and their reproduction in book form post- 
poned until theories are better established and more solidly supported by 
experimental results ? 

W. R. HAWTHORNE 


= 
¥ 
- 
| 


a 


JOURNAL OF 
FLUID MECHANICS 


Orders originating 


in 
UNITED STATES OF AMERICA 


and 


CANADA 


should be sent to the sole agents 


ACADEMIC PREss INC. 
111 Fifth Avenue, 
New York, 3, 
N.Y., U.S.A. 


M 
g 
¥ 
ts 
4 
4 
| 


JOURNAL OF FLUID MECHANICS 


Volume 2, Part 1 January 1957 


CONTENTS 


Dynamics of a dissociating gas. Part I Equilibrium flow 
by M. J. LIGHTHILL 


Transition processes in shock wave interactions 
by ROBERT G. JAHN. . 


The diffusion of smoke from a continuous elevated point- 
source into a turbulent atmosphere 


by F. B. SMITH 


Boundary layers whose streamlines are closed 
by W. W. WOOD 


A general solution of the equations of hydrodynamics 
by L. M. MILNE-THOMSON. . 


A boundary layer theorem, with applications to rotating 
cylinders 
by M. B. GLAUERT . 


Reviews . 


88 


7/ 


